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Introduction 


The scientific goal of the experiment is to test the equality of gravitational and inertial 
mass (i.e., to test the Principle of Equivalence) by measuring the independence of the rate 
of fall of bodies from their compositions. The measurement is accomplished by 
measuring the relative displacement (or equivalently acceleration) of two falling bodies 
of different materials which are the proof masses of a differential accelerometer. The 
experiment is a null experiment in which a result different from zero will indicate a 
violation of the Equivalence Principle. The final aim is to measure the Eotvos ratio 6g/g 
(differential acceleration/common acceleration) with a targeted accuracy that is much 
improved with respect to the state of the art (presently at several parts in 10 13 ). 

The estimate of the accuracy attainable with our experimental method is being carefully 
evaluated based on more refined analyses and advancements in the definition of the 
system. As the work progresses, we are setting our accuracy goal at improving the 
present EP test accuracy by about two orders of magnitude. 

The analyses carried out have focused on: (1) the evaluation of modal frequencies of the 
detector; (2) the dynamics of the instrument package after release following an imperfect 
release; (3) the computation of the higher-order mass moments on the proof masses 
generated by nearby masses; (4) thermal analysis of the instrument package during flight 
operation; (5) a preliminary design of the leveling and release mechanisms; and (6) a 
conceptual design of the cryostat and its interface. We are also cooperating with a 
cryostat designer on the topic of cryostat design development and costing. 

The project also involves an international cooperation with the Institute of Space Physics 
(IFSI/CNR) in Rome, Italy. The group at IFSI is in charge of prototyping the differential 
accelerometer and carrying out precursor laboratory measurements. Our Italian 
colleagues have received their external funding on this project from the Italian Space 
Agency (ASI) in November 2002. Thanks to some activity carried out in 2002 with IFSI 
internal funds and the new activities made possible by the ASI funds, they have been 
remarkably productive. They built a differential accelerometer laboratory prototype and 
carried out very important laboratory measurements focused on evaluating the common- 
mode rejection factor and on techniques for damping the natural dynamics of the 
accelerometer after release. The highlights of these activities are documented in a section 
of this report. 
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Update on Scientific Issues 


Need for Tests with Improved Accuracy 

The accuracy of the Equivalence Principle (EP) tests with laboratory proof masses on the 
ground is limited by the Earth's seismic noise and the weakness of suitable signal sources. 
Previous experiments include the famous torsion balance experiments of Eotvos (1890- 
1922)' " as well as the classical tests of the Equivalence Principle by Roll-Krotkov-Dicke 
(1964)'" utilizing a torsion balance which had an accuracy of 3 parts in 10" and the 1. 1. 
Shapiro et al. (1976) IV and Williams et al. (1976)' lunar laser ranging experiment with an 
accuracy of 3 parts in 10 12 . The present state of the art is at several parts in 10 12 both for 
lunar laser ranging tests (Williams et al., 1996) VI and for torsion balance tests (Adelberger 
etal., 1999)™. 

In two recent papers'" 1 lx , Damour et al., provide a theoretical justification, based on the 
string theory for expecting violations of the Equivalence Principle at an accuracy above 
one part in 10 12 which is slightly below the present state of the art in EP tests. Their 
argument centers on the inflationary rate of expansion of the early universe and on the 
decoupling of scalar fields from matter during the expansion of the universe. When the 
recent estimates of the early expansion rates are factored into the theory, they conclude 
that violations of the EP may well occur at accuracy much lower than previously thought. 

The conclusion of the analysis by Damour et al. gives strength to EP experiments with 
improved accuracy even if the improvement is not as dramatic as those promised by the 
space-based tests (of several orders of magnitude) of the EP. Another implication is that 
our goal should not be to push the accuracy as high as possible (and consequently 
increase the difficulty and cost of the experiment) but rather achieve a significant 
improvement in accuracy with respect to the present state of the art. 

One of the strength of our experiment is that by conducting it in free fall, the strength of 
the signal is increased by about three folds with respect to the signal available to 
experiments conducted in the laboratory. This simple fact simplifies the detector 
requirements vs. an equivalent-accuracy experiment conducted on the ground. However, 
in a drop experiment like ours, the measurement time is limited and, consequently, the 
detector must satisfy requirements associated with the ability to abate the transient 
dynamics after release. 

A rapid damping of the detector natural oscillations excited by the instrument release 
during the first few seconds of release is required. This point translates into the ability to 
switch the detector’s quality factor from the high value required during the data taking to 
low values during the damping phase. Our colleague at IFSI in Rome have already 
demonstrated with laboratory tests of a precursor accelerometer that this is indeed 
possible (see later on). 


7 



The test of the Equivalence Principle requires a differential measurement of acceleration. 
This fact has a positive consequence in terms of the rejection of accelerations that affect 
the two proof masses equally (common-mode type) and their effects on the differential 
acceleration. The capsule shields the free falling detector from external perturbations and 
provides a quiet environment for the experiment. However, the noise produced by the 
detector moving inside the co-moving capsule is not equal to zero. The residual gas in 
the evacuated capsule produces a tiny drag force and also allows the propagation of 
strongly-attenuated perturbations from the capsule’s wall. The evaluation carried out in 
the previous Annual Reporf on the noise related to residual gas in the capsule points to 
the need for a common-mode rejection factor of 10" 4 in order to meet the desired target 
accuracy. Our Italian colleagues have already carried out laboratory measurements on 
the instrument prototype and achieved the required common-mode rejection factor in the 
frequency band of interest. 
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Updated Reference Configuration 


Capsule mass and Size 

One of the important results of the analysis carried out is that the small vacuum 
chamber at the top of the capsule has been eliminated in favor of a fully cryogenic 
vacuum chamber (see Figure 1 ). The new solution eliminates the problem of the gravity 
gradient produced by the small cryostat in the proximity of the detector and also provides 
more clearance to the instrument package during the early stage of free fall. 




Spin and release 
mechanism 


Instrument package 
before release 


Cryostat 


Figure 1 Schematic of capsule with instrument package attached before release 
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The new configuration also has the advantage of using a more conventional design of 
the cryostat rather than one with doors opening at the bottom as in the old reference 
design. 

The preferred size of the vacuum chamber/cryostat was identified(see the 
analysis/optimization conducted in Annual Report#l) as a cylinder with internal diameter 
of 1 m and an internal height of 2.1 m which results in an external diameter for the 
capsule of about 1.4 m and a free-to-chose capsule length of roughly 5.6 m. Preliminary 
data on large cryostats (from Janis Research), indicate that the mass of a cryostat of the 
size considered above will fit well within the mass limit of the system. Such 
chamber/cryostat will allow free fall times in the range 24-28 s depending on the amount 
of ballast added to the capsule. 

The overall capsule mass could start at a minimum value of about a 1000 kg which 
provides a free fall time of about 24 s. The free-fall time can be increased by 1 s for 
every 250-kg ballast (or other equipment) added to the capsule up to the point where the a 
(limit) mass of 2000 kg is reached. The low-speed ballistic coefficient |3 () of the capsule 
will vary from about 6500 kg/m 2 for a capsule mass = 1000 kg to 13,000 kg/m 2 for a 
capsule mass of 2000 kg. Correspondingly, the free fall time will vary from 24 s to 28 s. 



Experiment Sequence 

Figure 3 shows the preliminary timetable of the experiment sequence. The 
experiment starts with the loading of the sensor into the vacuum chamber/cryostat about 
2 weeks before the planned launch. This operation is then followed by the pumping down 
of the chamber and the refrigeration of the sensor. After connecting the capsule to the 
gondola and the balloon, the balloon is launched. The estimated time to reach altitude is 
of order 3 hours. Upon reaching altitude, the attitude of the capsule is stabilized by the 
leveling mechanism on the gondola, the sensor is spun up, and the dynamics of the 
system is analyzed. When the dynamics is within the acceptable bounds, the capsule is 
released from the gondola and the sensor is released from the top of the chamber/cryostat 
immediately afterwards. The science data is taken during the free-fall phase in which the 
sensor spans the length of the chamber. Shortly after the sensor has reached the bottom 
of the capsule, the blut (first stage of the deceleration system) is released and, when the 
speed has decreased below the required value, the parachute is deployed. 



Figure 2 Preliminary timeline of experiment 



Updates on Capsule Dynamics 


Capsule Drop Dynamics 

The free fall time of the instrument package inside the capsule was computed for several 
values of the low-speed ballistic coefficient in the previous Annual Report. Figure 4 
depicts the capsule drop distance and Mach number vs. drop time. 
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Figure 3 Capsule drop and Mach number vs. drop time 


Following up the computation of the preferable range of values for the capsule mass 
and size. Figure 5 shows the relative distance covered by the instrument package during 
free fall with respect to the capsule for a low-speed ballistic coefficient ranging from 
6500 kg/m 2 to 13500 kg/m 2 . 

These results clearly indicate that a capsule of relatively compact size can provide a 
free fall time between 24 s and 28 s with a capsule mass ranging from 1000 kg to 2000 
kg, respectively. The higher value of the mass is for the same capsule with ballast. 
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Figure 4 Relative vertical distance vs. free fall time for limit values of interest of 
ballistic coefficient. Dots mark the expected operating points for the light capsule (M = 
1000 kg) and ballasted capsule (M = 2000 kg). 



Time (s) 


Figure 5 Relative lateral displacement due to wind shear vs. free fall time. Dots mark the 
expected operating points for the light capsule (M = 1000 kg) and ballasted capsule (M = 

2000 kg). 
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An internal diameter of the capsule of 1 m was also computed from the optimization 
process. This diameter determines the tolerance to wind shear acting sideways on the 
capsule. Figure 6 shows the lateral displacement of the capsule with respect to the free- 
falling instrument package vs. free fall time for the limit values of interest for the low- 
speed ballistic coefficient. The value adopted for the wind shear of 0.005 s" 1 is equivalent 
to a vertical gradient of 10 knots per km. This value is twice as high as the maximum 
wind shear reported" for the Air Force balloon base at Holloman, New Mexico. The 
lateral displacements due to wind shear are relatively small in the parameter range of 
interest. Furthermore, if the balloon is launched during the periodically-occurring wind 
reversal times (in April-May and September-October) the vertical wind gradient is much 
smaller than the value adopted for the computations shown here. 
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Modal Analysis of Detector 


Introduction 

The following is the modal analysis of the differential accelerometer for the EP test. 
The sensor is composed by two masses of different materials constrained to an external 
case by means of elastic springs. The measurement is performed by monitoring the 
displacement of the two masses with respect to the external case through the 
measurement of capacitance variations caused by the displacement. 

This analysis evaluates the eigenfrequencies and eigenmodes of the detector and the 
sensitivity of the modes to variations in the detector parameters. 

The model adopts six degrees of freedom, that is, one rotation and one translation for 
each body. The system parameters considered for the sensitivity analysis are masses and 
moments of inertia of the bodies, stiffness of the springs and their ratio. 

The influence of these parameters is evaluated on both the separation of the natural 
frequencies and the vibrational modes, particularly for what concerns the differential and 
common modes of the detector. 


Detector Dynamics Model 

The detector consists of two sensing masses having their Centers of Mass (CM) 
coincident with the CM of the external case. In the current detector model, one mass has 
a dumbbell shape, while the other is a hollow cylinder. The shape of the flight detector 
sensing masses will likely be different from this early choice. Nevertheless, the shape of 
the sensing masses, however, does not affect the generality of this study. The two 
sensing masses are constrained to the case by means of elastic springs, having a pivot 
axis parallel to the longitudinal axis of symmetry of the three elements. Figure 7 shows a 
schematic of the instrument model. 


B 


x 


spin axis 


pivot axis 


Figure 6 Schematic of detector model. 
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For this analysis, we assume that the capsule containing the sensor is fixed to the 
Earth surface, so that only the sensor is falling into it. This assumption does not affect the 
generality of the modal analysis. The Z axis is the axis of fall and of the Earth’s gravity 
acceleration. 

The reference system Cf is a frame parallel to the geographic X, Y, Z axes and 
displaced by the vector r z (see Fig. 8). Another coordinate system D is attached to the 
capsule and having axes parallel to the local normal frame. At this stage of the analysis, 
reference frame D is assumed to coincide with the local normal frame and the capsule is 
assumed stationary with respect to the Earth’s surface. This assumption doesn’t cause any 
loss of generality for the modal analysis. The coordinate that describes the fall of the case 
along the z axis is r z (t), where r is the position vector that identifies the body C (external 
case) with respect to the capsule and, finally, z indicates the axis of fall. 



Figure 7 Local normal (geographic) and Cf frames 

The body coordinate systems x, y, z are centered at the body CM, and fixed with the 
body itself. Each system has the x axis along the longitudinal axis of symmetry, 
corresponding also to the spin axis while the y and z axes are radial axes rotating with the 
body. 

A rotational and translational degrees of freedom have been considered for each 
mass. More precisely, the x rotation about the CM of each mass, and the y translation 
referred to the xyz frame. Therefore each body can rotate about the x axis through its CM 
and can translate along the y axis of our reference frame. Figure 9 shows the possible 
movements of each body. 
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Figure 8 Degrees of freedom of each body 


External forces 

The motion of the sensing mass consists of a rotation 0(t) about the pivot axis (related 
to the torsional stiffnesses of the springs) and a translation y(t) that depends on the 
flexural stiffness of the springs. The flexural stiffnesses are typically much greater than 
the torsional ones. The distance d in the figure identifies the distance between the CM of 
the sensing mass and the pivot axis. 

The external forces are written in the Cf reference frame, and consist of elastic forces. 
The elastic forces are produced by the springs connecting the sensing masses to the 
external case, and consist of a torsional torque and a force. In this analysis we assume 
that the rotational displacements are small, so that the displacement of the point S along 
the y axis due to a rotation around the CM is simply equal to 8d [i. e. we assume that 
cos(0) = 0]. 

The differentia] displacement of the point S between body A and body C and the 
elastic force along the y axis are therefore: 

Av a = y A - 0 A -d A - (y c - 0c'd A ) 

F yA = -k yA -Ay A = - k yA (y A - 0 A -d A - y c + 0 c -d A ) 

Where y A and y c , are the displacements along the y axis with respect to the Cf frame, 
0 A and 0 C are the rotations about the x axis of masses A and C, respectively. k vA is the 
flexural stiffness along y of the spring connecting body A to the case C. 
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The elastic torque about the CM of A is: 

N a = — k tA (0 A — 0c) — k yA -d A -Ay A = — k lA (0 A — 0 C ) — k yA - d A (y A — 0 A ’d A — yc + 0c'd A ) 
Where k tA is the torsional stiffness of the spring connecting A and C. 

Similarly, for the other two bodies we have: 

FyB = — k yB ‘Ay B = — ky B (yB — 0 b‘^b — yc ^c’^b) 

Nb = — k tB (0 B — 0c) ■ k yB ’d B -Ay B = -k tB (0 B — 0 C ) — k yB - d B (y B — 0 B *d B — yc + 0c'd B ) 

Fyc = — k yA (yc ■ 0 c -d A — y A + ^ A 'd A ) — k yB (y c — 0c‘d B — yB + ^b'^b) 

Nc = — k tA (0 C — 0 A ) — k yA -d A (y c - 0 c - d c — y A + 0 A 'd A ) — k tB (0 B — 0 C ) - k vB 'd B (y B — 0 B - d B 
— y c + 0 c 'd B ) 


These equations will be used in the next section to compute the stiffness matrix. 

Modal Analysis 

The mass and stiffness matrices are computed in this section. It has to be noted that 
out of 18 degrees of freedom (DoF) associated with the translational and rotational 
motion of the three masses only six of them are describing the motion of the degrees of 
freedom of interest. Specifically, only the translations along the y axis and the rotations 
about the x axis have to be considered for each body because the other motions do not 
affect the measurement. The state vector involved is, therefore, as follows: 

X = [}’A }’B yc ®A Qb &C ] 7 ( 8 ) 

which contains the translations along y axis and the rotations about the principal x axis of 
each body. 

The mass matrix is a diagonal matrix in a system of rigid bodies as follows: 


tn A 0 

0 nifi 

0 0 

0 0 

0 0 

0 0 


0 0 

0 0 

me 0 

0 ! Ax.x 

0 0 

0 0 


0 0 

0 0 

0 0 

0 0 

1 Bxx 0 

0 l Cxx 


(9) 


The stiffness matrix, obtained from the force and torque equations by using the 
displacement method, is: 


18 


0 


0 


k yA 

0 

~ k yA 

-d A k y A 

0 

d A kyA 


kyB 

—kyB 

0 


-d B k yB 

d B k yB 


-kyA 

—d A kyA 

CC 

’ 

1 

0 

kyA + kyB 

d A k yA 

d A k yA 

-) 

ktA + d A k yA 

d B k yB 

0 

d A k yA ~d B k vB 

s 

ri^t: 

*s3 

1 

1 


d A k xA 

d B k yB 

~d A kyA ~ cl skyB 
-k 


-d B ky B 
d B ky B 
0 

k rB + d B k yB 

0 ■) ■) 

~ k tB - ct Bk vB k tA + k tB + d A k v4 + d B k vB 


tA 


■d A k yA 

2 

-k t B ~ d B k yB 


( 10 ) 

We can observe that for an external body that is not constrained in its free fall, two 
rigid motions are possible, that is, the rigid body rotation around the x axis of the Cf 
frame, and the rigid body translation along the y axis. In order to carry out an analysis of 
the natural frequencies and natural modes of the system we can eliminate these two rigid 
motions so as to preserve only the elastic modes. This can be done by expressing one 
rotational and one translational coordinate of a body as a function of the corresponding 
coordinates of the other two bodies. 

A vector having the three translations equal to a constant and the three rotations equal 
to zero is a solution of the eigenproblem, so that this vector is an eigenvector. 
Consequently, all the others eigenvectors must be orthogonal to this one to yield: 

(<D 0 } t [M]{d>} = 0 

where <t> 0 is a rigid eigenvector, and <J> is another eigenvector. After evaluating the 
matrix products we obtain: 

OofmAO, + m B <I>2+ m c <t>,) = 0 

where now <t>j is the i-th component of a general eigenvector. Because O, = fx j5 where Xj 
is the i-th component of the vector x and f is a constant equal for all the components, we 
can also write 


m A y A + m B y B + m c y c = 0 

which allows to express one translational coordinate of a body as function of the 
translations of the other two. It can also be noted that in deriving the former equation we 
have obtained the conservation of the linear momentum for any elastic motion. 

We can obtain similar expressions for the rotational case as follows (from now on we 
drop the subscript “xx” in the moments of inertia): 

Ia 0 a + I B 9b + Ic 0c = 0 

We can now express the coordinates of a body as function of the analogous 
coordinates of the other two bodies. We can write therefore 

{x} = [C]{x c } 
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where x is the coordinate vector, C is a matrix, and x e is the new coordinate vector having 
the two coordinates of a body expressed as function of the others. In our case, we have 
eliminated the coordinates of body C, obtaining the following matrix [CJ: 


1 0 

0 1 

m A ™ B 

m c m c 

0 0 

0 0 

0 0 


0 0 

0 0 

0 0 

1 0 

0 1 

Ll 

! C J C 


( 11 ) 


and the vector X e being 
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( 12 ) 


The stiffness and mass matrices can be modified using the new coordinates by 
[KJ = [C] t [K] [C] 

[MJ = [C] t [M] [C] 

where M e and K e are the transformed mass and stiffness matrices which are not shown 
here for the sake of brevity. The dynamical matrix [D] = [M e ] '[K e ] that is needed to 
compute the eigenfrequencies of the elastic modes is as follows: 


1 


20 


b 

' 1 + t ' 


k yA m B 

d A k yA {l a + k : ) 

d A k yA I g 

k yA\ 

, m A ">C j 


in A nig 

I c m A 

I c m A 


kyB m A 

b 

— + — 1 

d B kyBl A 

dgkygd B + I c ) 


mgn\y 

k yB 

v nl C ) 

Ic m B 

Ic m B 

d A k 

yA ( m A + nl C ) 

d A k yA m B 

(Ct +I c){ktA + d A k y A ) 

h(k,A + d A k y A ) 


1 a"X; 



l A l C 

Vc 


d B kyB m A 

dgk Y g( nig + me ) ^ A ( + ^B^b j 

{ I B + l c)\kiB + (l Bky 


l B ni C 


I B nl C 

h l c 

hfc 


( 13 ) 

Once solved the eigenvalue problem for the constrained system (i.e. with the previous 
dynamical matrix), the total eigenvector including also the two variables eliminated can 
be obtained simply premultiplying the eigenvector obtained by the matrix [CJ. 

Eigenvalues and Eigenvectors. 

Using the dynamical matrix evaluated before an analysis of the eigenvalues and 
eigenvectors have been carried out in order to evaluate the influence of parameters on the 
natural frequencies and modal vectors. 

The parameters are the following: 

M a = mass of body A 

M b = mass of body B 

M c = mass of body C 

I A = moment of inertia of body A about the spin axis (x axis) 

I B = moment of inertia of body B about the spin axis (x axis) 

I c = moment of inertia of body C about the spin axis (x axis) 
d A = distance between the CM and the constrained point of body A 
d [f = distance between the CM and the constrained point of body B 
k, A = torsional stiffness of spring connecting bodies A and C 
k tB = torsional stiffness of spring connecting bodies B and C 
k A = bending stiffness of spring connecting bodies A and C 
k B = bending stiffness of spring connecting bodies B and C 
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Numerical values have been assigned to these parameters to evaluate natural 
frequencies and eigenvectors. Subsequently, the values of key parameters have been 
changed while keeping the others constant. 

The values used for the first evaluation are as follows: 

m A = m B = 5 kg 

I A = I B = 0.009 kg-nr (which corresponds to an external radius of body A of 0.06 m 
and a mass m A ) 

I c = 1 kg-nr 

d A = d B = 0.06 m 

k tA = k tB = 9.593 Nm/rad (corresponding to a first torsional frequency for a single- 
degree-of-freedom system of 3 Hz) 

k yA = k yB = 18000 ft 2 N/m (corresponding to a first lateral frequency for a single- 
degree-of-freedom system of 30 Hz) 

Using these values we obtain the following values: 

f, = 2.99 Hz; f, = 3.29 Hz 

f 3 = 52.13 Hz; f 4 = 55.22 Hz 

and the associated eigenvectors: 

v, = {0.0427736, -0.0427736, 0.705812, -0.705812} 

v 2 = {0.0326524, 0.0326524, 0.706352, 0.706352} 

v, = {-0.0209932, 0.0209932, 0.706795, -0.706795} 

v 4 = {-0.0210127, -0.0210127, 0.706795, 0.706795} 

Remembering that the eigenvector components represent in order y A , y B , 0 A , 0 B , we 
can make a few considerations about the natural modes. 

The first vector is composed of a translation of the masses and a rotation about the 
CM. The sign is opposite for the two masses, so that each of them moves in opposition to 
the other while the signs of y A and 0 A are the same so that the pivot point tends to remain 
in its original position. In the first mode body C is at rest, the motion being auto- 
compensating. In summary, the first mode corresponds to a differential mode of the two 
test masses. 




Figure 9 First natural mode (differential) 

The second vibrational mode has a similar type of motion, the only difference being 
that the two masses move in the same direction. In this mode body C moves in the 
opposite direction of the masses in order to equilibrate the motion. 



Figure 10 Second natural mode (common) 

The third and fourth modes involve the lateral stiffness of the springs as the sign of 
the rotation and translation of each mass is now opposite, so that the displacement of the 
constrained point is now the sum of the y translation and of 0 d A term (i.e., y A - 0 A d A ). 
The y translation depends on the lateral stiffness of the springs that is much higher than 
the torsional stiffness. The first two modes are analyzed in detail in the following. 
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Parametric Analysis 
K x influence 

The following plots show the variation of the first and second natural frequencies as 
k y varies. In each plot the components of the eigenvector is also plotted. 
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Figure 11 Dependence of first natural frequency on k y 
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Figure 12 Dependence of second natural frequency on k y 
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In the previous plots we can notice that the natural frequencies increase with k y , 
having as asymptotes 2.99 and 3.29 Hz respectively. Increasing the value of k y therefore 
contributes to separating the two frequencies up to the value of 0.3 Hz. The following 
plot shows the influence of k y on the displacement of the constrained point. 
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Figure 13 Dependence on k y of translation of constrained point in first mode 
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Figure 14 Dependence on k v of translation of constrained point in second mode 

It can be seen that as k y increases, the displacement of the point S (constrained point) 
decreases, so that the differential displacement in mode 1 is reduced. Also for mode 2 the 
same behavior can be noticed below. 
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Figure 15 Dependence of first natural frequency on I c 
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Figure 16 Dependence of second natural frequency on I c 

The influence of I c is reflected in the frequency of the second mode; decreasing the 
value of I c the natural frequency increases, yielding also to a larger separation between 
the two first frequencies. It can be observed that the separation of the frequencies 
decreases as I c increases, having as asymptote the value of separation of 0.3 Hz, 
corresponding to the separation provided only by k y for the specified k y and m c . I c has no 
influence on the displacement of the pivot point. 
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m c influence 


The effect of m c on the frequencies considered is to increase the separation of the 
frequencies as m c decreases, similarly to the effect of I c . The plot below shows this effect 
versus the m c /m A ratio. 


nat freq (Hz’ 



Figure 17 Dependence of natural frequency separation on m c /m A 

Combining the effect of I c and m c , the separation of the frequencies increases from a 
value of 0.3 Hz when I c and m c are much greater than I A and m A , to a value of 2.17 Hz 
when the ratios I C /I A approaches unity. 

Relative variation off, m A , k yA , k, A 

The parameters k vA , I A , k tA , d A , m A are now varied between the two sensing masses so 
to test the influence of possible imperfections in the mechanical and physical 
characteristics of the sensing masses and springs. The ratios p A /p B , where “p” indicates 
one of the above mentioned parameters have been varied between ±5% of the nominal 
value. 

The main influence of these parameters is on the difference of displacement of the 
two sensing masses in the common mode (mode 2), resulting in a differential 
displacement that is detected. The plot below shows the differential displacement (y A - 
y B )/y A versus the p A /p B ratio. 
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Figure 18 Influence of bodies A and B parameter’s ratio on the differential displacement 

As it can be seen, the most influent parameter is the distance of the constrained point, 
and subsequently the mass, and inertia ratio. This kind of dependence is explained by the 
fact that d A changes the moment of inertia around the constrained axis with a squared 
dependence; m A influences both the translational and the rotational motion, while I A 
affects only the rotational DoF. 

Concluding Remarks 

The analysis performed allowed to describe the differential accelerometer’s modes 
and frequencies and their dependence on the inertial and mechanical parameters. The 
system has been modeled using 6 DF, the rotation around the x axis and the translation 
along a radial axis for each body. 

The first two frequencies and eigenvectors have been analyzed in detail because these 
modes of rotation around the pivot point affect the acceleration measurement the most. 
The separation of these two frequencies depends mostly on the mass and inertia moment 
of the C body, and slightly on the lateral stiffness of the springs. The lower C mass and 
inertia are, the more the first two frequencies are separated (if I C /I A = 1 = m c /m A the 
difference is 2.17 Hz for a value of k y = 18000jt , while it is only 0.3 Hz if I C /I A = 1 10 
and m c /m A = 6); increasing k y contributes to reducing the displacement of the constrained 
point both for mode 1 and mode 2. 

A relative variation of the parameters ratio between body A and B has then been 
performed, to evaluate the influence of possible imperfections of the inertial and 
mechanical characteristics of the testing masses. 
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The differential displacement of the CM of the sensing masses has been analyzed. 
The k lA /k tB ratio is the most influent, causing a Ay difference of 5% as the previous ratio 
varies by 0.7%. The k yA /k lB ratio has the same influence. The ratio m A /m B has roughly 
half the influence of the spring stiffness ratio, causing a 5% Ay as the ratio varies by 
1.3%. The inertia ratio I A /I B produces a 5% variation in the y differential displacement as 
it varies by 2% from the unit value. Finally the lateral stiffness ratio k yA /k yB has a 
negligible influence compared to the others, i.e., approximately two orders of magnitude 
less, reaching a 0.025% Ay as the ratio varies by 1%. 
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Dynamics Analysis of Instrument Package/Detector 


Introduction 

The following analysis is a study of the accelerometer and instrument package 
dynamics. Two different designs of sensor have been taken into consideration as follows: 
(1) a sensor with two sensing masses and a case (three bodies); and (2) a sensor with one 
sensing mass (made of two different materials) and a case. In both types the measurement 
is carried out by monitoring the differential displacements between masses, through the 
measurement of capacitance variations caused by those displacements. 

In this study we derive the differential equations of motion for each sensor and carry 
out numerical integrations of cases of interest, starting from different initial conditions of 
the positions of the centers of mass (CM) and initial rotational velocities of the bodies in 
order to evaluate the influence of those parameters on the differential output of the 
accelerometer. 

Each body is defined using all its six degrees of freedom (dof) and elastic forces and 
torques are applied to each of them at the attachment points of the constraining springs. 

Description of sensor configurations 

Three-body sensor. 

The sensor is composed by two sensing masses having the CM coincident with the 
CM of the external case. One mass has a dumbbell shape, while the other is a hollow 
cylinder. These two masses are constrained to the case by means of elastic springs, 
having a pivot axis that is parallel to the longitudinal axis of symmetry of the sensor. All 
the bodies are spun about the symmetry axis x. Figure 20 shows a schematic of the 3- 
body sensor. 
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Figure 19 Schematic of three-body sensor 
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The motion of the sensor is essentially a rotation of the sensing masses about the 
pivot axis. The translational motion of the sensing masses, which depends on the high 
lateral stiffness of the springs is much smaller than the rotational motion. A rotation of 
the mass implies a variation of the distance between the sensing mass and the external 
case that changes the output capacitance. 

Two-body sensor. 

The two bodies sensor consists of an external case (body A) and an internal sensing 
mass. The internal mass is actually formed by two equal masses of different materials 
which are rigidly attached so as to behave like a single rigid body. The internal mass is 
connected to the external mass by means of elastic springs which have a lateral stiffness 
much higher than the torsional one. All the bodies are spun about the symmetry axis x. 

I 



Figure 20 Schematic of two-body sensor 

A violation of the EP during the fall would cause a rotation about the pivot axis of the 
sensing mass with respect to the sensor case. The measurement is then carried out by 
detecting the differential displacement (and therefore the capacitance variation) of the 
central disk of body B with respect to the four fixed capacitor plates attached to body A 
(see Fig. 21). The capacitor plates are arranged in such a way that a rotation of B about 
the pivot axis z unbalances the capacitive measurement bridge and produce an output 
signal. 

Reference frames 

A body reference frame has been attached to each body. The attitude of each body is 
identified by Euler’s angles, following the sequence 1-2-3 (i.e. by building the rotation 
matrix using sequential rotations about x, y, and z axes). 
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The relationships between the different reference frames have been written by means 
of rotation matrices and translation vectors. All the matrices are described in details in the 
following. The R 0A matrix transforms the inertial frame (X,„ Y 0 , Z 0 ) into the body A 
reference frame (X A , Y A , Z A ); 



Figure 21 Inertial and body reference frames 

The body frame is attached to the geometrical center of the body which, in ideal 
conditions, coincides with the CM of the body. The position of the CM (of body A for 
example) in the A reference frame is defined by (x CMA , y CMA , z CMA }. The coordinates of 
the origin of the A frame (R A ) with respect to the inertial 0 frame are: {x A (t), y A (t), z A (t)}. 

The matrix yielding the transformation from system X 0 , Y 0 , Z 0 to X A , Y A , Z A is 
composed by the following sequence of transformations: 

- rotation 0 xA (t) around X axis; 

- rotation 0 vA (t) around Y’ axis (transformed Y () axis after 0 xA (t) rotation); 

- rotation 0 zA (t) around Z” axis (transformed Z 0 axis after 0 xA (t) and 0 yA (t) rotations); 

- translation R A projected along the three axes X A , Y A , Z A . • 

The rotation matrix [R^] and the translation vector {R A } from system 0 (inertial) to 
system A (body A) are the following: 
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[*o/t] 


c(d yA )c(6. a ) c(d :A )s(0 xA )s(d yA ) + c(0 xA )s(6 :A ) -c(0 xA )c(d :A )s(d yA ) + s( 0 xA ),v( 0 :A ) 
-c(d yA )s(6 :A ) c(0 xA )c(6 :A ) - s(0 xA )s(0 yA )s(d :A ) c(0 :A )x(0 xA ) + c(0 xA )s(d yA )s(d :A ) 

s(d yA ) -c(0 xA )s(0 xA ) c(0 xA )c(0 yA ) 


where the superscript identifies the relevant reference frame. The translation vector is the 
R a vector, expressed in the A frame. This vector can be computed as follows: 

R a = [Roa] {-x a -y A -z a } t 

Where the vector {-x A , -y A , -z A > is the opposite of the vector {R A } expressed in the 0 
frame. 


We have: 



-c(0 yA )c(6 :A )x a - \c((L a )s( 6 xA ).s( 0 xA ) + c(0 xA )s(0 :A )|v A + [c{0 x/ , )cW :A )s(d xA ) - s(0 xA )s(0 :A ) )z A ' 
c(0 yA )s(0. A ).x A - [c(0 xA )c(0. a ) - s(0 xA )s(0 yA )s(0 :A )]v A - \c(0- A )s(0 xA ) + c(0 xA )s(0 vA )s(0 :A )]c. A 
-s(0 xA )x A + c(0 yA )s(0 xA )y A - c{0 xA )c(0 yA )z A 


For the opposite transformation, which leads from the body A coordinate system to 
the inertial frame, we use the transpose of the rotation matrix [R A() ] = [R 0A ] T and the 
vector is simply the R A vector expressed in the inertial coordinates {x A , y A , z A }. 

We have 


[M= 


c{0 xA )c(0 yA ) -c(^vA )s(0 xA ) 

C(0 :A )S(0 XA ) + c(0 xA )s(0 xA )s(0 :A ) c(0 xA )C(0 :A ) - s(0 xA )s(0 yA )S(0 :A ) 
-c(0 xA )C(0 :A )s(0 xA ) + s(0 xA )S(0 :A ) c(0 :A )s(0 xA )s(0 yA ) + c(0 xA )s(0 :A ) 


s(0 yA ) 

-c(0 yA )s(0 :A ) 

c(0 yA )c( 0, a ) 




x A 





In summary, we have the following transformations: 
Inertial to body A: 

[RoJ-lr 0 } + {R a A > 


Body to inertial: 

[Rao]- {r"} + {Ra°} 

The same procedure is followed for the other bodies. 
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Degrees of Freedom and equations of motion 

To determine the equation of motion for the different bodies the transformation 
matrices have been written for each body. The body coordinate systems are centered in 
the body’s geometrical center, and fixed with the body itself. The position of the CM is 
described using three coordinates (x CMA , y CMA , z CMA for body A). Each system has the x 
axis along the longitudinal axis of symmetry, corresponding also to the spin axis, the y 
and z axes are radial axes rotating with the body. 

Each body has all the six degrees of freedom, i.e. the three translations of the body 
frame origin, and the three rotations according to the Euler’s angles described above. 

To derive the equations of motion the translation accelerations of the CM and elastic 
forces have been written in the inertial system. The rotational equations are then 
projected onto the body axes. 

Translational accelerations 

In the following we describe the procedure for deriving the equations of motion by 
taking the two-body sensor as an example. The procedure to derive the equations for the 
three-body sensor is exactly the same with the added complexity of a higher number of 
equations. 

The CM vector in each body’s reference frame is: 

Body A: {r A } = { 

X CMA YCMA Z CMa} 

Body B: {r B } = {x CMB y CMB 

z cmb) 

The translational acceleration expressed in the inertial reference is obtained using the 
formula: 

{ a A } = [^Ao] ’ { r A } + ^[^o] ' } + [^Ao] ' {'A } + {^A } 

Where 

{a A } is the acceleration of body A in the inertial frame; 

[R ao ] is the rotation matrix from A to inertial frame; 

{r A } is the coordinates vector of body A CM in the A frame; 

{Ra°} is the translation vector expressed in inertial coordinates 

(°) indicates the time derivative and the dot (.) matrix multiplication. 

The same equations are used for the other bodies. 


34 



Elastic Forces 


To evaluate the elastic forces we define first the points where the springs are attached 
to each body. The position of these points are expressed in body frame. The distance 
between connected points is then computed from the actual position during the motion, 
projected in the A frame, and multiplied by the stiffness vector {k xA , k yA , k zA }. By 
following this procedure, it is possible to assign different stiffness to each degree of 
freedom. Finally the elastic forces are projected in the inertial coordinate system, and 
added to the equation of motion. The points connected by the springs are as follows (all 
distances are in meters): 

(Pai> = {0, 0, -75x1 0 y 

(Pa2> = {0,0, 75x1 O' 3 } 7 

{p B |} = {0, 0, -65x1 0'’} T 

{p I!2 } = {0, 0, 65x1 O '}'' 

The initial lengths of the two springs constraining the motion are as follows: 

{!.} = {Pa, - Pb,} = {0, 0, -10x10 y T 

{l2} = {p A 2-PB2> = {0,0, 10xl0-y 

The figure below shows the position of the points where the forces are applied 



(F | A } = [RJ • ([K] . ({p A1 } - [T 0A ].{T B0 .p B1 } - {1,})) 
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{F 2A }= [R ao ] • ([K] . ({p A2 } - [T 0A ].{T B0 .p B2 } - {1 2 })) 
{F A }={F IA } + {F 2A } 


In the above notation the operator “T ()A .x” indicates the combination of the 
multiplication for the rotation matrix [R<, A ] applied to the vector {x} plus the translation 
of the vector {R A A } which yields: 

{T oa .x} = [R oa ].{x} + {R a a } 

The expression T 0A .(T B0 .p B1 ) projects the coordinates of the point p B , (expressed in the 
body B frame) onto the body A coordinate system. Subsequently, the force is projected 
onto the inertial frame using the [R A0 ] matrix. A similar procedure is used to evaluate the 
forces acting on any other body pairs. 

The equations of translational motion for the two bodies A and B are: 


m A {a A }-{F A } = 0 


m B {a B }-{F B } = 0 


where m A and m B are the masses of body A and B, a A and a B are the accelerations, and F A 
and F b are the elastic forces: 


{F A } = {F IA } + {F 2A } 


{FJ = {F ib } + {F 2B }. 

Rotational accelerations: 

The angular velocity of each body is computed by using the rotation matrices that 
transform the coordinate system from the inertial to the body frame and conversely for 
the opposite transformation. The rotational velocity matrices of bodies A and B are 
derived by using the Cartan’s formula as follows: 

[ w a] = [KmH^o] 

[ w b] = [^ob]' [^bo] 

where [R 0A ] and [R A0 ] are the rotation matrices from coordinate system 0 to A and vice- 
versa. [R ob ] and [R b0 ] are the correspondent matrices for body B. 

It should be noted that the former expression lead to the skew symmetric matrix of 
the angular velocity from which the components of the angular velocity vector {co A } can 
be readily extracted. 

The rotational inertia forces are obtained using the formulas: 
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[ 1 b]' { c0 B } + [ C0 B ] ' [ l B ] ' } 


Where {co A } and [to A ] represent the angular velocity vector and matrix, respectively, 
for body A (and similarly for body B), [I A ] and [I B J are the inertia matrices for body A 
and B which, assuming principal axes, have the form 
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Elastic torques 

The elastic torques acting on each body are computed by using the definition of the 
attachment points and the expression of the elastic forces previously defined. The 
expressions of the torques in the body reference frame are the following: 


{T a } = {p A1 } x ([R oa ] . {F, a }) + {p A2 } X ([R oa ] . {F 2A }) 


{T b } = {p BI } X ([Rob] .{F 1b }) + {p H2 } X «R ob ] . {F 2B }) 


where x indicates the external product of vectors. 

The expressions of the rotational equations of motion are finally: 

+ [ w b]' \Jb\ { (0 b} ~{^b} = 0 


() 


The equations of motion for the three bodies sensor are evaluated at the same way as 
for the two bodies. The same quantities are derived and the same structure for the 
equations is followed. The following figure shows the positions of the constrained points 
for the three bodies sensor. 


The points p CAI , p CA2 , p CB , and p CB2 are the points of the external case C connected 
with the points p A1 , p A2 , p Bi , p B2 , respectively. In this case four springs are used to connect 
the two bodies to the case (two springs per each body). 
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Figure 23 Locations of attachment points 

The coordinates of the points are the following (each point is expressed in its own 
body system): 

{p A| } = {-40x1 O' 3 , -60x10 \ 0} T 
{p BI } = {-250xl0'\ -60x1 0" 3 , 0} T 
{PcAi} = {-50x1 O' 3 , -60x1 O' 3 , 0} T 
(Pcb i } = {-270x1 O' 3 , -60x1 0' 3 , 0} T 

{Ui) = {Pai — Pca I } 

Obi) = {Pbi - Pcb i } 

{p A| } = {40x1 0' 3 , -60x10 \ 0} T 
{p BI } = {250xl0'\ -60x1 0' 3 , 0} T 
{Pca i } = {50x1 0' 3 , -60x1 0 3 , 0} T 
{PcBil = {270x1 O' 3 , -60x1 O' 3 , 0} T 

{^a:} = {Pa2 — Pcaz} 

{^ 2 } = {Pb2 - PCB 2 } 
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Two Bodies 


The integration of the equations of motion has been performed using different values for 
the CM position and initial angular velocity of body A. 

A frequency analysis has been performed for each result, to evaluate the frequency 
content of the differential rotation around the z axis, which leads to the detection of a 
differential signal by the capacitors. 

The parameters in the equations are the following: 

m A = mass of body A 

m B = mass of body B 

I Ax = moment of inertia of body A about the x axis 

I Av = moment of inertia of body A about the y axis 

I Az = moment of inertia of body A about the z axis 

I Bx = moment of inertia of body B about the x axis 

I Bv = moment of inertia of body B about the y axis 

I Bz = moment of inertia of body B about the z axis 

k x = lateral stiffness of spring connecting bodies A and B (in A coordinates) 

k y = lateral stiffness of spring connecting bodies A and B (in A coordinates) 

k z = lateral stiffness of spring connecting bodies A and B (in A coordinates) 

k Hz = torsional stiffness of spring connecting bodies A and B 

x cma = x position of CM of body A (in body A coordinates) 

ycMA = y position of CM of body A (in body A coordinates) 

z cma = z position of CM of body A (in body A coordinates) 

x cmb = x position of CM of body B (in body B coordinates) 

ycMB = y position of CM of body B (in body B coordinates) 

z cmb = z position of CM of body B (in body B coordinates) 

The variables for the equations (dof) are: 

x A (t) = x coordinate of the geometrical center of body A (in inertial coordinates) 

y A (t) = y coordinate of the geometrical center of body A (in inertial coordinates) 

z A (t) = z coordinate of the geometrical center of body A (in inertial coordinates) 

x B (t) = x coordinate of the geometrical center of body B (in inertial coordinates) 

y B (t) = y coordinate of the geometrical center of body B (in inertial coordinates) 

z B (t) = z coordinate of the geometrical center of body B (in inertial coordinates) 

0 xA (t) = x rotation of body A 
0 yA (t) = y rotation of body A 
0 zA (t) = z rotation of body A 
0 xB (t) = x rotation of body A 
0 yB (t) = y rotation of body A 
0 zB (t) = z rotation of body A 
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Numerical values have been assigned to the parameters and then a study of the influence 
of the position of the CM with respect to the geometrical center of each body (CM 
parameters), also combined with an initial angular velocity of body B. The values used 
for the fixed parameters are in the following. 

We analyze in the following representative cases with prograde precession (i.e., ICx > 
ICy). The reason for favoring prograde precession over retrograde precession is that 
prograde precession brings about significantly smaller accelerations due to imperfect 
release (see later on) than retrograde precession. 

Two-Body Dynamics Simulation 

Prograde precession 


mA = 30 kg, mB = 1 kg, 

lAx = 95/100 kgm 2 , lAy = 146/100 kgm 2 , lAz = 146/100 kgm 2 , 
IBx = 1 7/1 00, I By = 1 7/1 00, IBz = 1 7/1 00, 
kx = 45000 N/m, ky -> 35000 N/m, kz -> 35000 N/m, 
kx2 = 45000 N/m, ky2 = 35000 N/m, kz2 = 35000 N/m, 
k0x = 1000 Nm/rad, k0y = 1000 Nm/rad, k0z = 54 Nm/rad; 

spin frequency -> 0.3 Hz 

initial vel theta z body A = 10' 2 rad/s 

initial nutation angle = 0 


xCMA = 0, yCMA = 0, zCMA = 0, 
xCMB = 10' 6 m, yCMB = 0, zCMB = 10' 6 m 


2. 99716 

7. 38688 

8. 08957 
42. 8045 
42.8045 
48. 5357 j 


Figure 24 Natural frequencies 
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Figure 31 Zoom of 0 zB - 0 zA FFT vs. frequency (Hz) 
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Three Body Dynamics Simulations 

The parameters in the equations are the following: 

m A = mass of body A 

m B = mass of body B 

m c = mass of body C 

I Ax = moment of inertia of body A about the x axis 

I Av = moment of inertia of body A about the y axis 

I Az = moment of inertia of body A about the z axis 

I Bx = moment of inertia of body B about the x axis 

I Bv = moment of inertia of body B about the y axis 

I Bz = moment of inertia of body B about the z axis 

I Cx = moment of inertia of body C about the x axis 

I Cv = moment of inertia of body C about the y axis 

I Cz = moment of inertia of body C about the z axis 

k Ax = lateral stiffness of spring connecting bodies A and C (in C coordinates) 

k Av = lateral stiffness of spring connecting bodies A and C (in C coordinates) 

k Az = lateral stiffness of spring connecting bodies A and C (in C coordinates) 

k Bx = lateral stiffness of spring connecting bodies B and C (in C coordinates) 

k Bv = lateral stiffness of spring connecting bodies B and C (in C coordinates) 

k Bz = lateral stiffness of spring connecting bodies B and C (in C coordinates) 

k HxA = torsional stiffness of spring connecting bodies A and C 

k HxB = torsional stiffness of spring connecting bodies B and C 

x cma = x position of CM of body A (in body A coordinates) 

Ycma = y position of CM of body A (in body A coordinates) 

z cma = z position of CM of body A (in body A coordinates) 

x cmb = x position of CM of body B (in body B coordinates) 

Ycmb = y position of CM of body B (in body B coordinates) 

z CMB z position of CM of body B (in body B coordinates) 

x cmc x position of CM of body C (in body C coordinates) 

y CM c y position of CM of body C (in body C coordinates) 

z cmc z position of CM of body C (in body C coordinates) 

The variables for the equations (dof) are: 

x A (t) = x coordinate of the geometrical center of body A (in inertial coordinates) 

y A (t) = y coordinate of the geometrical center of body A (in inertial coordinates) 

z A (t) = z coordinate of the geometrical center of body A (in inertial coordinates) 

x B (t) = x coordinate of the geometrical center of body B (in inertial coordinates) 

y B (t) = y coordinate of the geometrical center of body B (in inertial coordinates) 

z B (t) = z coordinate of the geometrical center of body B (in inertial coordinates) 

x c (t) = x coordinate of the geometrical center of body C (in inertial coordinates) 

y c (t) = y coordinate of the geometrical center of body C (in inertial coordinates) 

z c (t) = z coordinate of the geometrical center of body C (in inertial coordinates) 

0 xA (t) = x rotation of body A 
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0 yA (t) = y rotation of body A 
0 M (t) = z rotation of body A 
0 xB (t) = x rotation of body A 
0 vB (t) = y rotation of body A 
0 zB (t) = z rotation of body A 
0 zC (t) = z rotation of body C 
0 xC (t) = x rotation of body C 
0 yC (t) = y rotation of body C 


Prograde precession 

ICx > ICy 
spin freq = 0.3 Hz 


values: 

mA = 1 kg, mB = 1 kg, mC = 30 kg, 

lAx = 17/100 kgm 2 , lAy = 17/100 kgm 2 , lAz = 17/100 kgm 2 , 

IBx = 17/100 kgm 2 , IBy = 17/100 kgm 2 , IBz = 17/100 kgm 2 , 

ICx = 95/100 kgm 2 , ICy = 146/100 kgm 2 , ICz = 146/100 kgm 2 , 
kAlx = 45000 N/m, kAly = 35000 N/m, kAlz = 35000 N/m, 
kA2x = 45000 N/m, kA2y = 35000 N/m, kA2z = 35000 N/m, 
kA0x = 61 .68 Nm/rad, kA0y = 1000 Nm/rad, kA0z = 1000 Nm/rad, 
kBIx = 45000 N/m, kBIy = 35000 N/m, kBIz = 35000 N/m, 
kB2x = 45000 N/m, kB2y = 35000 N/m, kB2z = 35000 N/m, 
kB0x = 6168/100 Nm/rad, kB0y = 1000 Nm/rad, kB0z = 1000 Nm/rad 

initial vel theta z body C = 10' 2 rad/s 
Initial nutation angle = 0 

xCMA = 10‘ 6 m, yCMA = 0, zCMA = 10' 6 m, 
xCMB = 0, yCMB = 0, zCMB = 0, 
xCMC = 0, yCMC = 0, zCMC = 0 
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Figure 32 Natural frequencies 


The simulations have been run for 100 s to provide a good frequency 
resolution for the FFT plots. 



Figure 33 z A (m) projected in body C frame vs. time (s) 
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Figure 34 z A - z c (m) projected in body C frame vs. time (s) 



Figure 35 FFT of z A - z c vs. frequency (Hz) 



Figure 36 Zoom of z A - z c FFT vs. frequency (Hz) 
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Concluding Remarks 


The analysis carried out allowed to evaluate the influence of CM position and initial 
angular velocity errors in the dynamic of the 2 bodies and the 3 bodies sensor. The two 
systems have been modeled using 12 DoF, and 18 DoF respectively (six dof for each 
body). 

2- body results 

The position of the CM for body A (external case) did not affect the rotation about the z 
axis (measurement axis), while the initial nonzero condition on the angular velocity of the 
same body caused a differential signal to be detected. 

Analyzing the frequency content of the signal, though, no frequency peak was noticed at 
the spin frequency (0.5 Hz), while the natural frequency for the rotation about the z axis 
(3 Hz) and the precession frequency (0.18 Hz) were present. 

These frequencies can easily be changed modifying the torsional stiffness about z axis 
and changing the inertia ratio (IAz - IAx)/IAz, to which the precession frequency is 
directly related. 

It has to be noted that the spin frequency appears in the single 0zA or 0zB starting from 
an initial angular velocity of body A about y or z axes, but it is rejected when the 
difference between the two signals is performed. 

Also combining the initial conditions together (both CM positions and angular velocities) 
leaded to analogous results, resulting in a rejection of the spin frequency in the 
differential rotation. 

3 - body results 

Also in the case of the three bodies the position of the CM of the external case (body C) 
along the spin axis does not influence any of the variables; the displacement of the CM 
along the axes y and z cause the excitation of the y and z variables only, while the 
rotations are not affected. Moreover the differential displacement is not excited. 

An initial condition on the angular velocity (both along the y or z directions) excites all 
the dof of the system, both global and differential; the frequency pattern of the variables 
doesn’t show evidence of peaks at the spin frequency, neither in the global, nor in the 
differential component of the displacements. 

Applying all the different conditions together again all the dof are involved, and the 
frequency analysis evidences that the measurement differential signal (along z direction) 
doesn’t show any peak at the spin frequency, even if these peaks are present in other 
global variables (y and z) due to initial displacement of the CM of body C along y and/or 
z directions. 
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Detector Requirements Development 
Model for the Gravitational Perturbations Acting on Test Bodies 

Overview 


The test body (proof mass), falling inside the capsule, is subjected to non-negligible 
gravitational attraction by the capsule. In particular, the higher order gravitational 
potential plays a major role because of the elimination of the zero-order potential due to 
the equivalence principle. The model consists of a spinning test body inside a capsule: 
hollow cylinder covered with flat caps. The test body is released at the axis of the 
cylinder, and deviates from the axis during its fall. Therefore, our task is to compute the 
gravitational force and torque acting at the neighborhood of the fall. The fact that both the 
test mass and the capsule are closed finite bodies, increases the complexity of the 
problem. We took advantage of the size of the test mass; smaller than the cylinder radius, 
to derive an asymptotic analytical solution. We mostly concern with the force/torque at 
the modulation frequency. In other words, our model should evaluate the force/torque in 
a rotating body frame, rather than the static field in the capsule frame of reference. For 
that purpose we built a semi-analytical model. The main advantage of our model is its 
robustness. It can handle any configuration of test mass as well as any additional mass 
distribution in the test chamber. 

The purpose of this work is two folds. We need a working gravitational model for future 
simulation for the dynamics of the test body. We also need to choose the best inertia 
properties of the test body, in order to minimize the gravitational disturbance. 

The following report presents the analytical approach and the analytical results, and 
shows numerical investigations and the analysis of the results. In addition, we include 
few appendixes. Appendix A explains the computer code. Appendix B and C present 
analytical solutions for the gravitational attraction due to the cylinder and the caps, and 
check the validation of our numerical model. Appendix D proposes the approximation of 
the gravitational field and its gradient by using radial basis functions. 

General Approach 

The purpose of this work is to model and to analyze the gravitational attraction 
between the capsule and the test body. There are at least three ways to approach this 
problem. The first approach is to compute the force between each capsule mass element 
and a test body mass element, and to perform a double summation on these forces. If N is 
the number of capsule mass elements, and N B is the number of test body mass elements, 
then the cardinality of the computation is N ® N B . (The complexity of the computation 
is the product of the complexity of a single element computation by the cardinality). The 
second approach is a double integration over the bodies. The drawback of the first 
approach is the heavy computational effort, especially when we need an online 
computation in a dynamical process. The weakness of the second approach is in the 
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cumbersome analytical computation. Moreover, the integration will result in an 
asymptotic series of complicated functions. This will require intense computation. In 
addition, asymptotic series introduce a truncation error into the computation. We adopted 
a third approach for this model. 

Our computational approach is to consider the capsule as a discrete ensemble of 
lumped mass, where the resultant force and torque acting on the test body are the 
cumulative force and torque due to each capsule mass element. In other words, the 
interaction is between a finite body and a point mass. This approach is a tradeoff between 
the previous approaches, its cardinality is N . The main advantage is due to the flexibility 
of modeling any capsule shape, according to future requirements. 

The purpose of the following computation is to come up with simple closed-form 
analytical expressions, for the force and torque acting on the spinning test mass due to the 
capsule gravitational attraction. 

Gravitational Model 


The gravitational potential for finite size bodies, is: 


dM 

V = -G f f —dM 

9t r 
o B 


Where in our case, is the test mass (proof mass), and 9t 0 represents the attracting 
bodies, as the capsule (cylinder and caps) and the Earth. 

For the purpose of simplicity we will proceed with a representation of the test body, 
and a single element of the attracting mass, M ; - . 


The gravitational potential at a representative element mass M, due to the test body 

dM g _ 

is: V (M . ) = -GM . f , where r is the radius vector between an element mass 

1 X r 

of the test body and dM a as shown in Figure 53. 
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Figure 42 Gravitation model for test mass and outer attracting bodies 

Assuming a sufficient discretization, the total force acting on the test body due to all 

_ N 

outer mass elements is: F B a W(M,), where N is the number of outer mass elements. 

H 

The model utilizes two coordinate systems. The first is the capsule frame, denoted by: 
{ X,Y,Z } . X,Y are so far arbitrary, while Z coincides with the symmetry axis of the 
main cylinder. The test body frame is denoted by {x, v,z} , is attached to the test mass. So 
far, the origin and the orientation of the coordinate system are arbitrary. The gravitational 
potential will be expressed in body frame. The inertia coefficients of the body are 
constants in this frame. 

The analytical procedure is to express r as r = R- p , and to expand the potential by 

1 1 00 

the following power series: — = — S' 

r R n = 0 



n 

P„ 

( R-p\ 

UJ 

n 

R p J 


Here P n is the Legendre polynomial of degree n. 



The derivation consists of the following steps. Substituting the inverse radius 
approximation into the potential, carrying out the integration, while expressing vectors in 

body frame, as follows. R = (l e + me + ne_)R and p = xe t + ye t + ze , where 

{e x ,e ,e. } are body unit directions, and {l,m,n\ are direction cosines between R and the 

body axes. 


The resulting potential is a summation over all mass elements M s : 


V- l t 
/ = 1 


N GM 2 r "I 

^ ' [ / x + m y + n z J 


R 
GM 


R 


2 R 
GM 
2 R l 


i 2 2 2 

2~[( 3/ -l)J vv +(3m -1)J +(3n - 1) J ,, + 6(lm J YV + / n J xz + mil J v , )] 


p[/(5/“ - 3) J xxx + m(5m 2 - 3) J vvv + n(5n 2 - 3)J. 


+ 3 m (51 - 1) J YVV + 3// (51 - 1) J ry _ + 3/ (5m - 1) J YVV 


xx: 


-VVV 


2 2 2 
+ 3«(5 m - 1) J YVZ + 3/ (5» - 1) J xzz + 3m (5» -1)J yzz 


1 

+ 6 Imn J ] + O ( — t~) J 

R 5 


where j p q ,. are the inertia integrals defined as: 


j P a., = Sx p y q z r dm 
• ' m b 

For the second order, it is common to work with the moment of inertia. The following 
relations relate it to the inertia integrals: 


55 


% 


* 


J.X.X = 2^vv + A: -I A.v) 
Jyy = — ( I v v + I ZZ — l yy) 


*Av\ I. XX 

J.XZ = Arc 

Jyz = Ac 


We define 6 = {x, y, z} as the offset between the origin and the center of mass. It 
should be emphasized at this point that the origin does not necessarily coincide with the 
center of mass. Although the test bodies are typically axisymmetric, and ideally the 
geometrical center is the center of mass, a mass perturbation can divert the center of mass 
from the geometrical center. 

Note that if the origin is at the center of mass, the body frame coincides with the 
principal axes, and 0(1/ R A ) is neglected, then the resulting potential is the so-called 
MacCullagh formula: 




gm 1 m 2 GM 


j(trl-31fl), where l R =e R -\-e R = l\ xx + m\y y +n\ zz 

K 2 R 

is the projection of the second order inertia tensor I on R . 

The force acting on the body is the gradient of the potential (the negative gradient is 
the force acting on Mj. 


The resulting total force components, in body frame, are: 


4 
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F x = G 1 M ,{ M -4 + [(3 1, 2 - D.v + 3 /-in - v + 3 /- n - z 

i = 1 /?,• V 

+ ^4 L(f r-/ " I .vv “I.-: + 21,,)/,- + 3m,. (5/, 2 - l)I xy + 3/7, -(5/, 2 - 1)I, ; + 15/,-m, :U f l v; 
R; ^ 


+ ^3 [ '(35 A 4 - 30 /, 2 + 3 ) J, v , + f ( 7/»,. 2 - 3 )/,. m,. J vvv + f ( 7 n ,. 2 - 3)/,- », i ;:: 

Rj - L 

+ y ( 7/, 2 - 3 ) /,/!/,• J.vat + y( 7/, 2 - 3 )/,- /!,- J xv - + |( 35 / / 2 m / 2 - 5/ ; 2 - 5 m , 2 + D J AVV 
+ y ( 7 m ,- 2 - 1)/,- /?,• Jyy: + 2 ( 35/ ( 2 ;;, 2 - 5/, 2 - 5//, 2 + 1)J,,. + y ( 7//, 2 - 1)/,^ 


+ 3(7/,. - 1) m,-/? ,■ J.wc 1 + H.O.T} 


N Mr //?■ Mr r _ 0 _ _1 

F v = G I M, { y 2 + — f 3/,-m,- v + (3m,. - D.v + 3m,-/;,. ; 

i = 1 /?,■ /?,•■ L 


+ -^ [( : ^ r,. - I xx - I~ + 2I yy )///,. + 3/,. (5m, 2 - DI.w + 15/,-m.n,. I xz + 3 //,. ( 5 m , 2 - 1)1 v - 1 
2 

+ - y [|(7/, 2 -3)/, m,. J xxx + ^(35m, 4 - 30m, 2 + 3)J yyy + ^(7/7, 2 - 3)m,«,- J zzz 

+ |(35/, 2 m, 2 - 5m, 2 - 5/,- 2 + 1) J xxy + y (7/,“ - Dm,/?,J vv - + y (7m, 2 - 3) /,• m,. J xyy 
+ “(7m,. 2 - 3) «),/),. Jyy; + “(7»f - l)/,«i,J.v-; + |(35m, 2 /j, 2 - 5m, 2 - 5/?, 2 + D J yzz 


+ 3(7m- -I)/,-/?,- J xyz .] + A/.0.77 


„ A/ m b I M r 

F. = G 1 M,{— 1- + — f 

i = 1 */ V 


3/,- /?,• .v + 3m,/?,. v + 3(n,. - D z 


+ — ^4 [(fr,. - 1,, - I vv + 21.-)//,- + 1 5/,m,?,I xy + 3/, .(5//, 2 - W.,- + 3m,.(5//, 2 - 1)I V; ] 
R i 

+ — 3 -ly( 7/,- 2 - 3 ) /,- n,J vvv + |( 7 m , 2 - 3 )m,. //, J vvv + “( 35 / 7, 4 - 30 / 7, 2 + 3 ) J,.. 

R; ^ 


+ y (7/, 2 - l)/??,/i,J vw + “(35/,. 2 /7,- 2 - 5/, 2 - 5/7, 2 + 1) J„. + y (7m, 2 - D /,■ n,- J.wv 


+ -(35m,. 2 /?,. 2 - 5///, 2 -5/7, 2 + DJ VV - + y (7/7, 2 - D /,•/?,• J + y (7/7, 2 - 3) /71,/7,-J v __ 


15 ,, 2 


+ 3(7/?,. -DljUij J HOT) 
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where 


C 


A 


r. - ( I yy + I zz 2\ XX )L + (I xx + l zz 21 yy)m. + (I xx + I yy 21 zz ) n. 

Note that the force has the following order of magnitude: 

F ^G^^-{0{\) + oA + 0{^-) + ofe -) + H.O.T} 

i R- Ri Ri 2 R, 

Where L B is a typical length of the test body, « 1 , and — « 1 . 

Rj Rj 

The torque, acting on the test body, can be computed with respect to the center of 
mass or with respect to the geometrical center. Both coincide when the body has perfect 
symmetries. In reality there are mass imperfections that cause the center of mass to 
deviate from the geometrical center. 

Since we have no a-priori knowledge of the imperfections, our reference point for 
computing the torque is the geometrical center. The effect of the imperfection in 
demonstrated for a simple case. Let <5 be the radius vector between the reference point to 
C.M. The bias is due to a disturbing mass m P at a radius vector 8 P relative to the 
reference point. Therefore, the value of <5 is: 


m P 

mT 


Pp 


Thus, the effect of imperfections is proportional to its disturbing mass and to its 
deviation from the center of mass. 
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Figure 43 A Model for the Torque 

There are few approaches for computing the torque. Our approach is to integrate the 
torque due to mass element in each body. The torque with respect to an arbitrary point is, 
therefore: 

T () - Tc.m + <5 x F r 
where. 


TcM=Gf f dx^dM B dM 0 = J -T-dM B 

9t o&b r '=' r 

F b =G f f 4 dM B dM 0 = G / ^dM B 

m 0 d\ B r /=i m B r 

There are three approaches to compute the integrals. The choice of the approach 
depends on the size of each body. If both and are small compared with R, than 
we can expand r, and obtain an analytical approximation. If only is small, as in our 
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case, we can approximate 9t B but have to integrate numerically in di 0 . If both are large 
relative to r, then the only way is to perform a double numerical integration. In our case, 
not only is large, it is also surrounded by so there is no topologically way to 
expand r in 9^ . 

TV i 2 2 

T = 3 G 2 A/ | [ m n (/ -I ) + l n I - / m I + (n - m )/ ] 

x j _ i / ^3 i i zz yy i i xy i i xz i i yz 

i 

+ — ^-[« (5m 2 - 1)J - m (5n 2 - \)J + n (5/ 2 - 1)J +10 ImnJ 

2R i i vyy i > > i ny / / / xyy 

- in (5/^ - 1 )J — 10/ m n J + m (1 0 /?^ _ Sm^ + 1 )J 

/ / xxz i i i xzz i i i yyz 

- n (10m 2 - 5n 2 + \)J +10 l(n 2 -m 2 )J ] + U.O.T) 

i i i yzz i i i xyz 

N i 2 2 

T = 3G 2 M{ [/ n (I -I ) - m n I +lml + (/ - m )/ ] 

v / = 1 ‘ ' 1 xx zz i ' w / / yz i i xz 

i 

+ — (5/ 2 - 1)7 +l(5n 2 -\)J -\0lmnJ -n(5m 2 -\)J 

2 R 4 i i xxx i i zzz i i i xxv i i xvv 

-/ (10« 2 -5/ 2 + \)J + n (10/ 2 - 5n 2 + \)J +l(5m 2 -\)J +10 ImnJ 

i i i xxz i i i xzz i i yyz i i i vzz 

+ 10 m(l 2 -n 2 )J ] + H .O.T] 

i i i xyz 

N i 2 2 

T = 3G 2 M { [/ m (I -I ) - (m -1)1 +m n J -l n I ] 

z j _ i / ^3 / / yy xx i i xy i i xz i i yz 

i 

+ — — [m (5/ 2 - l)J -l{5m 2 -\)J + / (10m 2 - 5/ 2 + 1)J 

2R 4 i i xxx i i vvv i i i XXV 

I ... 

- m (1 0/ 2 - 5 m 2 + 1 )J + 10/ m n J + m (5/7 2 - 1)J - 10/ m n J 

i i i xyy i i i xxz i i xzz i i i yyz 

- I (5/7 2 - \)J +10 n(m 2 -l 2 )J } + n.O.T} 

i i yzz i i i xyz 

The order of magnitude of the torque is: 
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T * {CK-^tt) + 0{^\) + H.O.T) 

i K, K- 

If the origin is at the center of mass and the body frame coincides with the principal 
axes, then the resulting torque is reduced the well-known expression (based on 
MacCullagh formula): 

N M 

7>3G(I 22 -I„)2-> : h, 

/=l A, 

N M 

Ty «3G(I„-I 22 ) 2-fli',! 

i= I A, 

N M 

T : .3C(I vy -l„)2-J-/,m, 

/ = 1A," 

We would like to propose a more general approach for a closed-form solution of the 
force/torque in terms of their frequencies. 

Let us expand the force in potential orders as well as a Fourier series in the spin 
angle, assuming a pure spin around the axial axis of the test body. 

F = Gf -r-j 2 {«o + j, [a { n d ] cos inp + b^ sin ncp]}dM 0 
:it t ,R d = o /i=i 

where d is the degree of the inertia. Let J ,• = / x' y J z k dM B than d = i + j + k. 

Assuming that we have carried out the integration (analytically or by numerical means), 
F can be expressed as: 


_ ^max d+\ _ _ , , , 

F = G 2 {^q + 2 \-K? ] cosiup+ B ( ,l } sinnq)]} 

< 1=0 n-1 

;;(</> u(</) 

where = / " r dM , and 5 (c/) = f —— dM . 

91„R‘'* 2 

These coefficients may be viewed as a generalization of the inertia coefficients of9i 0 . 
If 9 ? 0 is topologically connected (i.e., is outside 9f 0 ), and if is sufficiently 
small, then one may asymptotically expand these integrals. This will result in the inertia 
coefficients of9i 0 . In our case, these coefficients can be computed numerically, or 
sometimes analytically, when 9f 0 is simple enough. 

For a more general rotation, one needs to express the direction cosines in terms of the 
general transformation matrix (rather than the single-axis rotation). 


61 



Order of Magnitude and Similarity Analysis 






The dominant sources of the gravitational attraction are the capsule, the Earth, and to 
a lesser degree the Moon. Let us perform a qualitative comparison of each source. First, 
we should note that the force acting on the test mass is not a homogeneous function of the 
mass and the radius (as a force acting on a point mass). This is because the parameters 


F oco( Mb \ + 0\ 

f m b ) 

{ F 3 j 


Mb \ 


(M b ) 

\ + 0\ 

+ 0 

\R 2 ) ' 

1 ! 

\ R 4 ) 


{ R 5 J 


, or 


F a F (0) + F (l> + F^ 2) + F (3) , where F (0) represents the 2-Body term, and so on (the 
order enumeration is according with the corresponding Inertia order). Therefore, it is 
impossible to express the ratio of forces from different sources as 

F(M\,R\)/F(M2,R2)<xf((M\/Ml) m ,(Rl/Rl)"). 

On the other hand, each order is homogeneous, that is 
F (k) {M\,R\)/ F (k \M2,R2) = (Ml/M2)(R2/Rlf k) 

Therefore, in order to compare different sources of attraction, we must consider each 
order separately. 

The following qualitative discussion concerns the force difference between the two 
test bodies. Let us denotes /()as a generic function, then a generic expression for the 

force acting on test mass B. is: 


\R‘ 


F 1 =f(Mf,)0 — + f(p.;)a— + /(/ 2/ ;) 


\R 


R q 


+ f(hp)0\ 



We present two models for the force difference. In the first model the two test bodies 
are centered, and the inertia components of the bodies are different. In the second model, 
the bodies are of equal inertia, but are not centered. These two models are idealization of 
the two sources of difference, while the reality is probably a combination of these. 

The force difference for the first model is: 


AF = F 8 2 -F B ' = /(AM;) 


R 2 / 


+ /(Ap;) 



+ /(A/ 2 ;)0|— +/(A/ 3 ;) 

R 1 



The force difference for the second model is represented by a differential. For central 
deviation e , the force differential is: 
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dF*F B i - F B 1 s VF-f=/(M;)o(^] + /(p;)4^) + /(/ 2 ;)o(-^U/(/,;)d 

\r- J U / U' / 


( 1 ' 
V/? 6 , 


The ratio between force terms of order k from two different sources 1 and 2, is: 


F' k '(M ] ,R l ) M, | 

(r 2 \ 

F' k) (M 2 ,R 2 ) ~ M 2 

v / 


The following table compares the force order of magnitude due to different sources. 
Note that we eliminated order 0 (2-Body), since the test body is in free fall. 


Table 1 Ratio of force terms for few orders for the dominant gravitational sources 


Inertia order 

Capsule / Earth 

Capsule / Moon 

Earth / Moon 


4 *10’’ 

8 * 1 0 4 

2 * 10 7 

2 

3 * 10 J 

3 * 10" 

l*10 l) 

3 

2*10" 

1*10 22 

6*10"’ 

4 

1*10" 

4*10 50 

4* 10 I: 


It is obvious that the effect of the Moon is negligible. The Earth is dominant only 
with respect to the gravity gradient. Again, as the order under consideration becomes 
higher, the effect of the capsule increases. 

To gain more insight into the nature of gravitational attraction let us examine the 
force. It can be rewritten in the following form: 

F * °(f) o(c)+ 0 (f) 0<?< °(f) 1 o(e " ) + °(fH‘ ) 

Where c represents the direction cosines. The notations c° and c L show the nature of 
the degree of the direction cosines (odd or even). For example, a square of a particular 
direction cosine has an even degree, while a product of the three direction cosines has an 
odd degree. Now, suppose that the test mass is at the center of the cylinder. Since for 
each direction cosine to a mass element there is an opposite direction cosine, all the odd 
terms are cancelled out. Therefore, the only non-zero contribution is due to the even 
terms. 

Another important engineering issue is the role of the test body mass and size. For 
this purpose we perform a similarity analysis, as shown in the following. 
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We formulated the gravitational potential as an asymptotic series in \/R n . Each term 
in the series has an inertial tensor of degree n - 1. For example, the first term is of order 
n = 1 and its inertia tensor is of order 0 (the test body mass). The third term is of degree 
n = 3 and its inertia tensor is of order 2. This non-uniformity introduces a problem, as 
mentioned earlier, when one computes the forces for a particular test mass, and wants to 
deduct the forces for a different scale of test mass. Mathematically speaking, if F(M^) 

and F(Mfi 2 ) are the forces due to two different test masses, then the ratio between the 

111 

forces is not a homogeneous function, that is, F(M ) / F(M B2 ) * / M B2 ) ) 

where g is a function and m is the degree of the homogeneity. Note however that each 
term of the series is homogeneous by itself. For example, the first term is homogeneous 
of degree one, that is, the ratio between forces equals the ratio between the masses. 

We will distinguish between two cases. In both, the mass distribution of the different 
test masses can be scaled. The first case is when the difference in masses is due to a 
different density. Since the forces are homogeneous of degree one in density, the ratio 
between the forces is equal to the ratio between the masses. In the second case the density 
is the same, and the different mass is due to different sizes. In this case we will apply a 
different similarity to each term. Let L and M be scaling factors for the length and the 

1/3 

mass, respectively. If the density of different test bodies is the same, then L^M . Let 

m = n - 1 be the order of the moment of inertia. Its similarity dimension is L"'M, or 

M l+m/ \ or L"' +3 . The acceleration similarity is M"' /3 or L m . The following table 
summarizes the similarity dimension for each term in the potential. 
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Table 2 Similarity relations for the gravitational potential terms 




Inertia 

Similarity 

Acceleration 

Similarity 


0 

M (V) 

1 (1) 

2 

1 

M 4 ' 3 ( L 4 ) 

M ,n (. 1 ) 

3 

2 

m 5/3 (l 5 ) 

M 2 ' 3 (L 2 ) 

4 

3 

M 2 (L 6 ) 

M (I 3 ) 

5 

4 

M ln (L 7 ) 



To summarize, given the force (or the acceleration) due to a particular test body, we 
can deduct the force (acceleration) due to another geometrically scaled test mass. The 
acceleration ratio (of two different test masses) versus similarity dimensions (mass and 
length) is illustrated in the following plot. 
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Figure 44 Scaling of Acceleration of Test Mass 
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The maximum potential order considered in this analysis is the hexadecapole (n = 5) 
which involves the fourth-order inertia integrals. We computed numerically the forces 
associated with the hexadecapole for a body with equal second-order inertia integrals (in 
order to minimize the quadrupole term). The results show that, at the worst expected 
location inside the capsule, the maximum acceleration associated with the hexadecapole 
term for a 1-kg proof mass are less than 10' 16 g. In other words, for the accuracy goal of 
this experiment there is no need for belted cylinders (which reduce the hexadecapole 
component even further). Proof masses with equal second-order inertia integrals (or 
alternatively moments of inertia) and construction accuracy of order a few microns are 
sufficient to make the contributions of all the higher-order gravity terms negligible. 


* 

Effect of Inertia imperfection on the Spinning Test Body 


t 


p 


Our goal is to minimize the gravitational forces acting on the test mass. Observing the 
force equations reveals that the dominant term under our control is the term 
corresponding to the second order inertia. The offset effect may be bigger to start with 
but it can be reduced by centering. The direct way to minimize the forces is to require 
equal second moments of inertia. The left over forces is due to the imperfection, and to 
higher order inertia integrals. The purpose of the following discussion is to explore the 
effect of each term for a spinning test body. 


First we will evaluate the mixed-inertial terms due to the imperfection. 


We assume that the imperfection is due to disturbance in mass and in length. The 
nominal test body is a perfect cylinder, with radius R B and length L B . In order to minimize 
the force we require that I xx = Iyy = l X z- This constraint results in a given aspect ratio of 

the cylinder as follows: L B = t[ 3R b . Given the mass density, the mass or the sizes of the 


cylinder are now functions of a single free parameter. For example, given the mass, M B , 


and the density, p B , the length is: L B = 


3 Mr 


[jip B 


The similarity dimension of the k-order inertia is: [/] = Ml } . Thus, the perturbed k- 
ordered inertia is: [d/J = L k 6 M + kML k ~ ] 6L , or 
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6M' 

+ k 

6L 

— 

= 



/. 


_M . 


L 


Next we consider the effect of the spin. The test body is spinning about its x-axis with 
a frequency that will be regarded as a fundamental frequency, or PI (period one). We are 
primarily concerned with PI which is the period of the measured signal. For that purpose 
we will analyze/^ by substituting the direction cosines, shown below, into the force 
expression. 
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/ = / (t\P 1); m = m (t;P 1); n = n (?;P1) 
i i i i i i 

It is important to note that the forces are expressed in body frame. The direction 
cosines are PI -periodic. Moreover, if the body fall such that the body x-axis is parallel to 
the capsule X-axis, only m and n are periodic, as will be seen from the following 
discussion. 

The periodicity in time is introduced through the direction cosines that represent the 
orientation of the body frame relative to the capsule frame. The transformation between 
the two is: 



The resulting direction cosines are: 

/=- = — + 1 ( 0 ) 

R R 

v Ycd + ZsG 

m = — = 

R R 

z -YsQ + Zct) 

" ~ R ~ R 

I is a cyclic function only if the body deviates from the X-axis. If the test body is 
perfect, then only the point mass term and the first term of 0(\/R A ) contribute to Fy . 
Otherwise, m,n introduce higher harmonics. If the power of the direction cosines is even 
(N e ), than the additional harmonics are: . For odd power (N () ) the 

additional harmonics are: P ] ,P i ,"-,P fi . The dominant time dependency of a perfect 

body is introduced through the first order attraction on a point mass. It results in a 
fundamental frequency, because the gravitational attraction is static in capsule frame, 
while the measurement is in the rotating body frame. 

The expected frequencies are the following. PI from the point mass termO(l/7? ;: ). 
PI and P2 from the offset term o(l//?’ ). Since T. contributes P2, the term of o(l/7? I * * 4 ) 

contributes PI and P2. However, because of the almost axial-symmetricity of the test 
body, T. is almost a constant, so the dominant frequency of o(l//? 4 ) should be PI. The 

term of o(l/P' ) contributes PI, P2, P3 and P4, where P2 should be dominant. 

The analysis above is depends whether the x-axis is parallel to the X-axis. If it is, 
then / = constant . Otherwise / = 1(6) and higher frequencies would be involved. In this 
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case, higher frequencies would be invoked. The explicit orientation dependency of period 
one becomes: 


-Body 

r 

v 


(6;P.)= G 2 m\ 
i = 1 




+— 

R 3 - v 


+ J_p< 4 ) 

R 4 y 


+ -o <5) 

R i y 


1 


where 

F^"'* = M (Y cos 6+Z sin6) 
v B i i 

=3 M n X (Y cos d + Z sind)x 
y B ' i i 

F (4) =\-[AX 2 Y, +-(3 B + C)(Y 3 + Y.Z 2 )]-BY lcos0 
v [2 / / 4 i ii if 

+ \-[AX 2 Z + —(35 + C)(Y 2 Z . + Z. 3 )] - BZ . \ sin 6 
[2 t i 4 / i i if 

F (5) = -{IXf- - 3 )(Y cos 6 + Z, sin0)X/_ m . 
y 2 

and X i - XI /?, Y = Y / R, Z i = Z / R are the direction cosines in the capsule frame. 

An interesting case is the effect of the attraction of a point mass located on the Z- 
axis. It may represent the Earth attraction on horizontally falling body, or the caps on the 
top or the button of the capsule. 

Let the attracting mass M , be at Zl = +F, in the plan X = Y = 0 . Using the relations: 
Z, = R ] sgn(Z,), Z, = sgn(Z,), Z =1, and substituting the attracting mass coordinates 
results in the following simple equation: 


F^° dy (9;P v l = 0 )= GM,\ 


M R I 1 

— 7- + -(75 + 5C) 

R 2 8 R 4 


sgn(Zj )sin0 


Again, the relations above are the terms corresponding to PI. The other frequencies 
may be observed from analytical expansions, or from a numerical frequency analysis, as 
demonstrated in the next section. 
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Numerical Investigation of an Imperfect Spinning Test Body 


The scenario under investigation is shown in Figure 56. The capsule consists of a 
hollow cylinder, covered by two flat caps and a total mass of 500 kg. The test body mass 
is 1 kg with sizes: =0.04 12m and L B =0.07 13m. The geometrical and mass 

inhomogeneity errors are: 6L/L = 10" 4 and 6M/M = 3 O' 5 . Consequently, the offset 

components are: x = 1 .7 x 10~ 6 m , y = 2 x 10 _6 m , 5 = 1 .7 x 10~ 6 m. The components of the 

inertial matrix are: I xx = 8.5 x 10 -4 , I rr = 8.501 x 10 -4 , l zz = 8.5015 x 10~ 4 . Note that 
the nominal values of those are the same. The perturbations prevent the elimination of the 
following inequalities: / xx + fyy — ^ ^ xx + ^ zz — ** 0, /yy ^zz — 2/^y ^ 0. 

l XY = 3x IO -8 , I xz = 3 x 10 -8 , I YZ = 4 x 10~ 8 . The components of the third-order 
inertia tensor are: I xxx = 5 x 10 -10 , l YYY = 8 x I0~ 10 , l 777 = 8 x IQ" 10 . I xxy = 6 x I0 _l °, 
l xxz =6xlO“ 10 , l XYY = 7xl0‘ ! °,/ ryz =8xl0“ l() , I xzz =6xi0”'°, I YZZ = 8 x IO" 10 , 
I xyz = 7x10 10 . 

tz 


2.3 [m] 



1.2 [m] 



* 


t 


Figure 45 Coordinate systems for the test body and the capsule 
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The gravitational force and the torque acting on the test mass were computed at 9 
points inside the capsule. Three points along the axial axis, three points along an axis 
deviated by 0.1 m from the axial axis, and three points along an axis deviated by 0.2 m 
from the axial axis. The ninth point is the farthest away, thus representing the worst case. 



-1 - 0.8 - 0.6 - 0.4 - 0.2 0 0.2 0.4 0.6 0.8 1 


Figure 46 Tested Points inside the Capsule 

We are mostly interested in F y ( F y is in phase with F z ) and in T x . Because of the 

free fall, we exclude the 2-Body terms. The following plots show the angular history of 
the force/torque for a single revolution with respect to the x-axis, and the corresponding 
Fast Fourier Transform. The latter is normalized, such that the highest amplitude is one. 
The results are given in the following 18 plots. The analysis of the frequency spectrum 
« requires a closed look at the equations for the force and torque, considering the particular 

location of the test mass relative to the capsule. Some of the results are non intuitive. The 
4 general approach to estimate the results is to specify each term according to the degree of 

the direction cosines (odd or even), and consider possible cancellations due to anti- 
symmetricity with respect to the radius vector between the test mass and each capsule 
element. The results are shown in the following figures. 
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Figure 47 Ordered Force and harmonics at Point 1 due to Capsule Attraction 
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Figure 48 Ordered Torque and harmonics at Point 1 due to Capsule Attraction 
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Figure 49 Ordered Force and harmonics at Point 2 due to Capsule Attraction 
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Figure 50 Ordered Torque and harmonics at Point 2 due to Capsule Attraction 
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Figure 51 Ordered Force and harmonics at Point 3 due to Capsule Attraction 
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Figure 52 Ordered Torque and harmonics at Point 3 due to Capsule Attraction 
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Figure 53 Ordered Force and harmonics at Point 4 due to Capsule Attraction 
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Figure 54 Ordered Torque and harmonics at Point 4 due to Capsule Attraction 
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Figure 55 Ordered Force and harmonics at Point 5 due to Capsule Attraction 



spin angle [rad] frequency [HZ] 


Figure 56 Ordered Torque and harmonics at Point 5 due to Capsule Attraction 
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Figure 57 Ordered Force and harmonics at Point 6 due to Capsule Attraction 
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Figure 58 Ordered Torque and harmonics at Point 6 due to Capsule Attraction 
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Figure 59 Ordered Force and harmonics at Point 7 due to Capsule Attraction 
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Figure 60 Ordered Torque and harmonics at Point 7 due to Capsule Attraction 
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Figure 61 Ordered Force and harmonics at Point 8 due to Capsule Attraction 
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Figure 62 Ordered Torque and harmonics at Point 8 due to Capsule Attraction 
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spin angle (rad) frequency (HZ) 

Figure 63 Ordered Force and harmonics at Point 9 due to Capsule Attraction 
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Figure 64 Ordered Torque and harmonics at Point 9 due to Capsule Attraction 
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The frequency spectrum for the force and the torque are shown in the following 
tables. The dominant frequencies are in bolded face. 

Table 3 Frequency spectrum along the capsule axis 



First Order 

Second Order 

Third Order 

F 

y 

P0, P2 

P0, PI, P2, P3 

P2, P4 

T 

x 

P2 

P2 

P2, P3, P4 


Table 4 Frequency spectrum to the side of the capsule axis 



First Order 

Second Order 

Third Order 


P0, P2 

P0, PI, P2, P3 

P1,P2,P3,P4 

T 

X 

P2 

P2 

P1,P2,P3 


Table 5 Frequency spectrum close to the cap 



First Order 

Second Order 

Third Order 


F 

PI, P2 

PI, P2, P3 

PI, P2, P3,P4 


T 

X 

PI, P2 

P2 

PI, P2, P3 



The orders in the above tables are the following. Order 1- offset or inertia order 1, 
0(1 //? 3 ) in force. Order 2- inertia order 2, degenerates to MacCullagh formula for a 
perfect body (or with respect to the principal axes), 0(\/R 4 ) in force. Order3- inertia 
order 3, 0(1/ R 5 ) in force. 

It is important to emphasize that the relative orders depend on the mass imperfection. 
For example, we should not rush to a conclusion regarding the contribution of the offset. 
The dominant contribution for a perfectly symmetrical body is from the second order. In 
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this case the torque is in agreement with the MacCullagh formula. The spectrum related 
to the second order can be deducted by noticing that the direction cosines, m and n, have 
period one: m -» P \ , n -* P\ . Therefore, 

Tx = (Izz - Iyy ) m n -* P2 , Ty ~ ( Ixx - Izz) n —* P \ , Tz ~ ( lyy - lxx ) m -» PI . 

Note that Iyy~Izz, thus both the offset term and MacCullagh term are amplified by 
the offset. Tx is null for a perfectly symmetrical body. Ty and Tz show various 
frequencies. However, Ty and Tz acting on a perfectly symmetrical body contribute PI 
according to MacCullagh formula. Also note that order 4 may play a role as well, because 
we assumed equal principal inertia of order 2. This assumption nulls the perfect body 
contribution of the second order but it does not cancel completely the term of order 4. 

The frequency PI is dominant only at the extreme location (point 9), due to the offset. 
We investigated the vicinity of this point, and found that the amplitude of Pi decreases as 
we move farther from the cap. 


Concluding Remarks 

We presented an analytical and a numerical analysis of the gravitational perturbations 
acting on the test mass due to the capsule. The outcome of this work is a closed-form 
formulation for the force/torque as well as a robust and interactive computer code. 

Our analysis predicts that the gravitational perturbations, acting on the test mass due 
to the capsule attraction, are safely within the limit required by the experiment. 
Moreover, along most of the test mass trajectory, the perturbations are far smaller than 
the limit, and their frequencies are different from the modulation frequency. 

The conclusions from our investigations, concerning the test masses are: 

The test masses should be smaller than about 10 cm and lighter than about 2 kg. 

The second order principal moments of inertia should be equal within 
construction tolerances 6L/L < 10 4 . 

There is no need for belted cylinders (as used in the STEP satellite experiment) 
for the accuracy target of this experiment. 

- The test masses should be as perfect as possible in terms of their mass distribution 
(5M/M- 10 4 ). 

The investigation carried out is essential for the definition of the tolerable sizes, 
masses, characteristics of the moments of inertia, and construction accuracy of the 
sensing masses. 
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Release/Leveling Mechanism Development 


Introductory remarks 

The leveling and release mechanism must be able to release the detector with initial 
conditions that do not impair the accuracy of the science measurement. Once the effects 
of the higher-order mass moments on the proof masses are made negligible through the 
right shapes and sizes, the most important external perturbation is the Earth's gravity 
gradient. Other perturbations are of course present but they can be made negligible 
through thermal design, reducing the pressure inside the detector, and shielding the 
detector from magnetic disturbances. The noise components associated with the Earth’s 
gravity gradient manifest themselves at twice the spin frequency and (depending on the 
orientation of the spin axis) at the spin frequency. The latter are the damaging 
components. The gravity component of importance to us is the g xv where x is the spin 
axis and y the sensitive axis of the accelerometer. As shown in Ref. x, this component is 
proportional to the product <J)6 X where (j) is the elevation angle of the (body) symmetry 
axis with respect to the horizontal plane (defined by the local gravity) and 6 X is the 
centering error along the spin axis between the CMs of the proof masses. The formulas 
derived in Ref. x for the Earth’s gravity gradient components will be utilized at the end of 
this section to set a limit on the tolerable tilt angle at release. 

Additional harmonic components come from the rotational dynamics of the 
instrument package. These harmonics are related to the inertia characteristics of the 
package and the rotational velocity errors at release. In summary, the leveling and 
release mechanisms must provide an orientation of the spin axis close to horizontal (to 
reduce the Earth’s gravity gradient component) and rotational velocity errors sufficiently 
low to avoid saturation of the accelerometer output. 

The detector has also its own elastic dynamics (as shown in previous section) which 
is excited by the conditions at release. We can conservatively assume that the release 
will excite the elastic dynamics of the detector up to its end of scale. The experiment 
strategy is to damp the elastic oscillations through electrical dissipative forces (see later 
on for experimental results) exercised for a few seconds after release. Once the elastic 
(natural) oscillations are abated to a level well within the dynamic range of the 
instrument, the electrical dissipative forces are removed and the detector operates as a 
high-Q detector. After the natural oscillations are abated, the oscillations of the proof 
masses will be forced by the rotational dynamics of the detector during the fall. 

In order to understand the effects of initial errors at release upon the differential 
accelerometer output, we need to develop a simplified model of the accelerometer. This 
model must contain the key dynamical elements but must also have a sufficiently simple 
formulation that shows analytically the origin and frequency content of the proof masses 
acceleration. 



Simplified dynamical model 

The differential accelerometer consists of proof masses that are about twenty times 
lighter than the mass of the instrument package. Moreover, the proof masses are 
expected to move with respect to the CM of the whole package by only microns during 
the early transient phase and by nanometers during most of the fall. Consequently, the 
rotational dynamics of the package is unaffected (and this will be confirmed in follow-up 
analyses) by the tiny motions of the proof masses. On the contrary, the rotational 
dynamics of the package drives the motion of the proof masses. In summary and with a 
good degree of approximation, the rotational dynamics of the package can be assumed to 
be steady. The solution of the attitude motion of the rigid body is available in close form 
for a free spinning body. We need to right the equation of a proof mass that is 
mechanically constrained to the rotating (and wobbling) instrument package in order to 
understand the origin and frequency content of the acceleration output of a single proof 
mass. The analysis can be readily extended to two proof masses. 

The general expression for the acceleration of a test mass at point P with respect to O 
in a rotating system Fean be written in matrix form as: 
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where [oj] is the angular rate matrix of the reference system and {x y z} 1 the position of 
the point P with respect to the center of the reference system O (placed at the instrument 
package CM). 

For a rigid body the equations are projected onto the axes attached to F itself (body 
axes) where the position vector is assumed fixed with respect to the center of F. 

Consequently, eqn. (7) transforms into: 
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Where the angular rate matrix projected onto the body axes can be computed from the 
rotation matrix that relates the body axes to an inertial system: 
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( 3 ) 


After substituting eqn. (9) into (8), we derive the well-known acceleration matrixes 
that provide the acceleration gradient at the point P with respect to O as observed by the 
rotating observer (or equivalently a proof mass of the detector): 
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Euler acceleration matrix 


Consequently, the (apparent) acceleration gradient tensor (in body axes) between a 
point P and the system CM at O is: 

[A] = [C] + [E] (10) 

In order to compute the expression of to x , oo r w z we need to consider the dynamics of 
the instrument package that houses the detector. In other words, we must solve the Euler 
equations for the free-torque case of a rigid body. We will consider an inertially 
axisymmetric body and regard x as the longitudinal axis of inertial symmetry. The non- 
symmetric body also has a known solution but it is more complicated as it involves 
elliptic integrals. 

After setting v = component of the angular rate along the axis of symmetry, the well- 
known solution of the Euler equations is obtained as follows: 


tu v = v 

• co v = -co, cos(Qr) (11) 

to. = to, sin(Qt) 

in which the origin of t is at the time when (o y attains its maximum value and where: 

Q = — — — v (body precession rate) 

h 
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transverse angular velocity 


I 2 2 I 2 2 

(0,=^Jc0y +C0- =*\J(Oy Q +W- 0 

The quantities to y0 and to z0 are the initial components of the angular velocity 
orthogonal to the symmetry axis. In our case, they can be interpreted as the rotational 
velocity errors at release of the instrument package. 

After substituting eqn. (11) into the matrices [C] and [EJ, we obtain the component of 
the overall acceleration gradient matrix: 

a xx = -co t 

a xx = a xx = -o), cos(£2t)(£2 - v) 

a x - =a. x = -to, sin(Qf)(Q-v) 

a xx = -^-co 2 [l -cos(2Q/)] + v 2 

1 1 

a x: = a :x = -—oij sin(2Qr) 

a. : - ^eo 2 [l + cos(2Qr)] + v 2 


The acceleration vector a measured at P and projected onto the body axes, is simply: 

M-MM (i2) 

where {6} = [6 X , 5 y , 6J 1 is the position vector from O to P (we have changed the notation 
to the 6s to highlight the fact that we are dealing with very small distances). The 
component of the acceleration orthogonal to the symmetry axis (i.e., along the sensitive 
axis of the accelerometer) is as follows: 
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a. = -co t (v - Q)cos (Qt)d x + —to, sin(2 Qt)d x + 


+ —cof( 1 + cos(2 Qt)) 


6 . 


( 13 ) 


This acceleration is the dominant acceleration experienced by the proof mass once the 
natural oscillations have been abated. In fact the amplitude of the residual natural 
oscillations can be made orders of magnitude smaller than the magnitude of the 6 vector 
through the initial damping. Equation (13) highlights several important points as follows: 
(1) the rotational velocity errors at release (encapsulated into oo t ) combine with the 
displacement error 6 to produce an acceleration output along the sensitive axis; and (2) 
this acceleration component is modulated by the body precession rate and not by the spin 
frequency. This fact was observed previously in the results of the general elastic model 
of the detector. The conclusion is very important because it implies that the errors at 
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release do not impact directly the harmonic component of a possible EP signal violation 
(that appears at the spin frequency). The overall advantage of the gyroscopic body (with 
non-spherical inertia ellipsoid) is that the precession rate is different from the spin rate. 
Moreover, the precession rate can be chosen at will by selecting the ratio of the moment 
of inertia of the instrument package in order to be less damaging to the ability of 
extracting the signal from the noise. Since both co p and 2 co p appear in the expression of 
the acceleration, it is advisable to make 0 ) p non commensurate with v (the spin rate). The 
selection of the moments of inertia ratio (and hence co p ) and the ability to extract a signal 
with a strength at the threshold sensitivity of the detector will be carried out in next year 
analysis. 


* 


Derivation of requirements 


The first requirement derives directly from the analysis carried out in Ref. x and it is 
related to the strength of the components of the Earth’s gravity gradient at the spin 
frequency. The resulting acceleration component (from Ref. x) is 


da, < 3-^cos(</>)sin(</>)sin(vr)<5 v + — ^-cos 2 (0)[sin(2vr)(5 v + cos(2vf)<5 r j 
D " 2 R~ 


R 

JL 

R 3 


^COS 2 (0)-1 


(14) 


for an instrument package spinning about x, a sensitive axis along z, and <j) the elevation 
angle of the symmetry axis with respect to the horizontal plane defined by the local 
(gravity) vertical. Equation (14) highlights the three components produced by the Earth’s 
field at the spin frequency, at twice the spin frequency, and a dc term. The strength of the 
gravity gradient component at the spin frequency must be less than the threshold signal at 
SxlO' 13 g. A product (J) 6 X of 0.1 deg-micron will safely meet the previous condition. 
Consequently, we can either be more relaxed on the leveling/release mechanism (i.e., (|) < 
1 deg and 6 X < 0. 1 micron) or on the centering between the two proof masses along the 
spin axis (i.e., ()) < 0.1 deg and 6 X < 1 micron). Both options will be kept open for the 
time being because they involve several technical considerations related to the 
mechanization of the leveling/release mechanism and the calibration of the differential 
accelerometer. A choice between the two options will be searched for after a dynamic 
analysis of the release mechanism and laboratory tests on the differential accelerometer 
prototype to assess the difficulties involved with the accurate centering of the proof 
masses along the spin axis. 

The requirement on the rotational velocity errors at release can be readily obtained from 
the following considerations. First, the rotational velocity error must be sufficiently 
small not to saturate the instrument output. In this case we do not have to worry about 
the instrument sensitivity because with appropriate inertia characteristics there will be no 
components at the spin frequency related to the precession dynamics of the instrument 
package. In this case the stronger component is associated with the first term of eqn. (14). 
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A rotational velocity of 1 deg/s will only produce a signal of order 10 '' g which is well 
within the dynamic range for the expected values of spin and precession frequencies. 

However, we must also worry about the amplitude of the precession (i.e., the nutation 
angle) because through coning of the body axis the accelerometer will sense the Earth's 
gravity gradient. The frequency sensed through coning is not exactly the spin frequency 
because the coning adds a frequency modulation at the precession frequency. 
Nevertheless, we will assume conservatively that the amplitude of the nutation angle 
must be less than previously indicated for the gravity gradient related noise. For an 
axisymmetric body the amplitude of the body-axis coning (i.e., the nutation angle) is as 
follows: 


tan(0) = 1l^L (15) 

h, v 

where co, is the transverse angular velocity at release, I., and 1, are the moments of inertia 
about the transverse and the symmetry axis (x-axis in our case), respectively. If we 
assume an upper bound for the nutation angle of 1 deg, a spin frequency of 0.5 Hz and a 
(worst-case) inertia ratio I/I., ~ 3, we obtain co, < 1 deg/s. Likewise, if we had assumed a 
maximum 0 < 0. 1 deg then we would have obtained to, < 0. 1 deg/s. In summary, we will 
establish the following two sets of requirements 

a) 6 X < 1 //m (centering error between proof masses CMs along the spin axis) 

(j> < 0. 1 deg (verticality) 

co, < 0. 1 deg/s (angular rate error at release) 

b) 6 X < 0. 1 /cm (centering error between proof masses CMs along the spin axis) 

<() < 1 deg (verticality) 

co, < 0 deg/s (angular rate error at release) 


We will start by base lining the release and leveling mechanism according to option a 
while keeping option b open if the tight centering of the proof masses of the detector 
turns out to be feasible. 
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Mechanical Report 


Release sequence 

The release of the package takes place in 2 stages. Initially the overall instrument is 
released from the balloon. Then, very soon after, the sensor package is released, in two 
steps, within the vacuum chamber. 

Overall Instrument Release 

The combination of wind shear during ascent, the extremely low air resistance at float, 
and the limited time at float make it necessary to assume that the balloon will be rocking 
at the time of package release. Since it is required that the package be vertical, and non- 
spinning through out the experiment, the release strategy must account for this. The 
present concept, developed in the reporting period, is to actively separate the orientation 
dynamics of the instrument from the behavior of the balloon. That is, we have placed a 3- 
axis gimbals between the instrument release mechanism and the portion of the gondola 
that stays with the balloon. 

As the balloon rockets back and forth, and rotates, the gimbals will be driven to ensure 
that the instrument package angular orientation is unchanged. The sole effect being that 
the package will be translating back and forth, and moving up and down slightly. During 
this process the linear motion will be tracked from the ground. In order to ensure that we 
release the package when it is as close to unaccelerated as possible, we will release it at 
one of the extremes of the balloon pendulum swing. 

Since the sensor release does not have any orientation adjustment capability, the overall 
instrument release mechanism must be designed to meet all the alignment and stability 
requirements that cover the release of both the instrument capsule and the sensor. 

Sensor Release 

The sensor release concept involves a two stage support. The first stage is a simple spit 
that holds the sensor package through release of the overall instrument from the balloon. 
Once the capsule is in free-fall the sensor release mechanism pulls back far enough to 
disengage the spit, and draw it to a position that the sensor will not hit it on the way 
down. At this point the sensor is held between 6 springs, 3 on each side. These springs 
will be sized so that they can not impart more force than is allowed by the rotation 
stability requirements. Any transient sensor motion induced in the spit-support release 
will be allowed to damp out in the 6 spring. Once the instrument package is spinning at a 
stable rate, and oscillation from the initial stage of the release have damped, the 
mechanism second stage will pull back further, and the instrument will drop. 

The next stage of the work on this design will be to create a detailed simulation of the 
entire system, in order to examine its behavior in detail. 
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Einstein Elevator 
3-Axes Gimbals Mechanism 


Figure 65 Schematic of leveling mechanism 



Einstein Elevator 
Instrument Package. Spin-Up 
& Release Mechanism 


Figure 66 Instrument package spin-up and release mechanism 
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Einstein Elevator 
Instrument Package. Spin-Up & 
Release Mechanism 


Figure 67 Front view of release mechanism 

Stabilizing Spring (3) 



Release Spindle Engaged To 
Instrument Package 


Figure 68 Detail of release mechanism (with spindle engaged) 


90 






Stabilising Spring (3) 
In Free Length Mode 



Release Spindle Disengaged From 
Instrument Package 




a 


Figure 69 Detail of release mechanism with spindle disengaged and springs still engaged 


Thermal Design 

One of the key advances this year was in the area of the thermal design and increased 
detail in the thermal model. In order to assist in this area the main supporting structure for 
the sensor package had to be sized. Though most of the instrument design remains at the 
concept level, it became important to examine the design forces that controlled this 
component’s dimensions. The sizing involved balancing the need to minimize the heat 
transfer, and therefore the support’s cross sectional area, while at the same time minimize 
the support deflection, or therefore the support’s structure moment. Studies of viable 
cryogenic materials suggested that a Kevlar composite would be the best material to use. 

Since we intent to use a composite material, and by the nature of the design, one with a 
thin wall, a further constraint on the design was a limitation on the allowable 
compression load. Coupling these requirements we were able to size the main shaft in % 

such a way that is can support the instrument without allowing out-of-spec instrument 
deflection, while dropping thermal conductivity well below the allowable level. ♦ 

Below are the present support parameters: 

o Length = 0.12m 
o Diameter = 0.05m 
o Wall thickness = 0.00 1 2m 
o Kt/m = 1 ,9x10 s W*m/K, at 10°K 
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Cool down Techniques 

Though this will be covered in more detail in the thermal section, the preferred method 
for lowering the instrument temperature to 10K has been re-examined and changed. 
Initially we had examined using radiation alone to reduce the instrument operation 
temperature to 10K. It was quickly clear that this approach was inadequate. Next we 
examined having a cold strap in place prior to launching the balloon. We would then rely 
on just radiation from the cold Dewar walls to maintain the temperature. Though this is 
effective, it introduces issues related to manipulating the cold strap. 

To avoid this we have since base lined a third approach, one involving cooling the 
instrument down by flowing cooled GHe through the vacuum chamber until the 
instrument is cool. We would then pumping out the remaining GHe from the inside of the 
Dewar. This give us a fast way to cool the instrument without having to manipulate a 
cold strap. 

Dewar Operation Techniques 

One of the larger instrument dynamics stability issues that needs to be addressed is the 
effects of the liquid cryogens during free fall: 

o the fact that they may slush during the fall 

o the fact that they will continue to boil during the fall 

o the fact that they needs to be vented, which will cause a propulsion effect. 

In order to avoid these effects we are beginning to examine novel Dewar operation 
techniques. The standard approach for maintaining temperatures at this level is a Dewar 
with a vacuum shroud and 2 cryogen shrouds, LN2 and LHe. But the LN2 is heavy. It can 
slush, it tends to boil similarly to water, and when it is vented it produces some 
propulsion, all of which will cause some dynamic reaction in the Dewar that could 
disrupt the experiment. The LHe, on the other hand is quite light and is far less likely to 
cause a dynamic disturbance. 

For these reasons we have decided to examine operating the Dewar, at least through the 
experiment phase, without LN2. The first approach that we examined was allowing the 
LN2 to completely vent prior to dropping the experiment package from the balloon. This 
is an acceptable approach, though it does produce some limitation on system timing. We 
would rather be able to begin the experiment when the system is settled, rather than 
waiting longer for the LN2 to boil off. Instead we are going to examine the possibility of 
allowing the LN2 boil off on the ground, and then flowing the boiled off LHe into LN2 
shroud. 

Dewar Layout 

We have moved the Dewar opening from the top of the chamber to the bottom. The 
initial Dewar design had the opening for inserting the instrument package into the top of 
the Dewar. This meant that all the feedthroughs into the system, as well as all the service 
feedthroughs for the Dewar itself were broken anytime the Dewar was opened. In 
addition, the load path from the balloon to most of the instrument mass goes through this 
connection. By moving the opening to the bottom, all the feedthroughs remain intact 
when loading the system, and only the instrument crash protection is supported through 
the Dewar cover (bottom plate). 




Figure 70 Cryostat top detail with feedthrough interfaces 



Thermal Analysis 


The simplified model discussed previously has been further developed, particularly in the 
area of isolation of the sensor from the (relatively) warmer preamp and electronics 
module. In fact, the electronics has to be kept above a minimum temperature (for this 
model, we have assumed -65C) to prevent damage to the components, and possibly 
needs to be warmed further to a “cold-start” temperature at the initiation of experiment 
calibration. 

The design consists of a sensor supported by a Kevlar-reinforced shaft (discussed 
previously) that is suspended by the rotation fittings across the dewar diameter. We have 
assumed the electronics module and preamp are balanced disks supported by this same 
shaft, both to one side of the centrally located sensor module. Between these various 
elements we have places radiation shields to minimize views of the warmer elements to 
the sensor. The model representation is shown in Figure 71. 

The electronics are shown in red as a result of their maintained temperature; the preamp 
is just visible between the two largest disks. 



Figure 1 Thermal model representation 



We have assumed that during the entire process (including cool down of the experiment) 
that the electronics is maintained above -65C, as illustrated by the red disk at the left end 
of the figure. 

There are two fundamental conductive paths between the various elements in the model: 

the supporting shaft and the wires connecting the electronic elements. The hollow shaft 

has been sized to support the masses with sufficient stiffness to allow the experiment to 

be balanced. The wires (currently assumed to be 4 pairs between each module) have 

been sized to provide less than 1-ohm resistance per circuit. Pure-metal wire like copper 

or silver have sharply peaked thermal conductivity in this temperature regime (10-25K) 

as well as rapidly changing electrical conductivity, so we have assumed for this analysis 

that the wires are Constantan, which has very little change in electrical characteristics j. 

from room temperature down to 4K, but its thermal conductivity drops by a factor of 100, 

providing effective thermal isolation at the experiment temperature. The conductors in 

the model have been sized to accurately represent the shaft and wires, and all the 

important materials are modeled with temperature- varying conductivity and specific heat. 

It has been shown in earlier studies that the experiment cannot be cooled effectively by 
radiative exchange with the dewar alone - it requires the addition of a thermal strap to 
achieve the desired starting temperature in a reasonable period of time. However, a 
mechanical strap (or straps) adds several complications to the design. The largest mass is 
the sensor itself, yet this is the least desirable location to affix a strap. We are currently 
considering an alternative approach: to flow gaseous helium over the experiment to cool 
it. This has the advantage that after starting conditions are achieved; the helium can be 
pumped out, allowing us to take advantage of the effective isolation of the mechanical 
hardware. 

We have run a study of the warm-up of the preamp and experiment with nominal power 
inputs and the previously mentioned initial condition of the electronics module at 208K 
(-65C). Using an assumed cooling of the preamp based on area ratio to the sensor and 
heat flow through the wires, the equilibrium temperature of the preamp prior to startup is 
approximately 20K. Using these initial conditions, and wires and support tube modeled 
as described above, the temperature rise of the preamp and the temperatures and 
estimated gradients in the sensor are shown in Figures 72 and 73. 
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Figure 3 Temperature rise at various locations of sensor 
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This study has yielded several significant conclusions for the design of the experiment: 

• Properly sized alloy wires provide sufficient electrical conductance without 
significant impact on the thermal performance 

• The most significant thermal conductance is the shaft supporting the sensor and 
other modules. 

• Radiative heating of the sensor can be well controlled with concentric disks 
between modules. These disks, if of significant thickness, also provide thermal 
mass to slow the heating of the sensor. 

It should be noted that while the temperature rates of change fall within the previously 
described limits (current estimate of maximum temperature in sensor is 0.0004 K/sec), 
improvements can be contemplated which would control these rates better. For example, 
the model currently has the shaft connected directly to the center of the aluminum endcap 
of the sensor module, essentially a worst-case solution. A set of Kevlar straps between 
the shaft and the sensor could provide high stiffness with lower conductance. However, 
such improvements, as well as better fidelity modeling of the inner structure of the 
sensor, require more sensor design detail. 
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Figure 4 Sensor temperature distribution at end of 20 min. 
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Instrument Electronics overview 


The overall communication and control system for the experiment is shown below. 
Minimum electronics beyond that already in the detector are added to provide a data 
collection and status link to electronics mounted outside the vacuum chamber. The link 
electronics on the detector is a battery-powered RS-232 to IR transceiver comprised of a 
110 ° field of view IR emitter-detector pair driven by an IR encoder-decoder and linked 
to an identical pair on the side of the chamber. The total power required for this portion 
of the interface should be less than lOOmW. Presently the electronics will need to run at - 
40 ° C, but we will investigate if we can go lower. 

The power and telemetry control box will be a comprised of a Book-sized PC with an 
Ethernet interface to a procured telemetry radio, such as the Freewave FGR1115RE 
900MHz spread spectrum radio. A small IR to RS-232 module provides the interface to 
the IR transceiver inside the capsule. The externally mounted system will run from a 
single +12VDC 85 W power source. The unit is expected to occupy a volume of 
40cm*30cm*25cm. 

The gondola will contain a standard SIP module for ground communication and 
control. We will need to provide an interface between the SIP telemetry interface and a 
radio that links to the experiment. This interface will essentially be the same design as the 
power and telemetry control box in the instrument. 

Other means of communication within the capsule were considered, including a low 
power 2.4 GHz radio link, and this option could be a viable alternative if the IR link is 
not feasible. The major advantage to this approach would that the link does not depend on 
position, and the power would be about the same. 



Figure 5 Schematic of telemetry links 
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System Requirements Development 


Table 6 is an update of key requirements of major subsystems. 


Table 1 Development of Requirements vs. Design Drivers (version #2 


Design Driver 


Free-fall (time > 20 s) 


Amplifier noise (white) 


, ; - V. v , 

1 ' • . -'%■ . 

ir-v&rr? 



Brownian noise (white) 


Temperature gradients 
[Radiometer effect, (to)] 


Acceleration noise inside 
capsule in free fall 


Earth’s gravity gradients 


Cryostat’s gravity 
gradients 
(distributed mass) 


Gravity gradients of lump 
masses on board capsule 


Magnetic disturbances 


Higher-order mass 
moments 



AT/Ax < 0.1 K/m T uniformity 
inside cryostat 


p c < 10 mBar 


Structural and 
attitude 
freqs. » to 


< 10" 12 g (2(0) Centering of 

< 10 15 g (©) proof masses (8 X ) 


< 10 12 g (2©) 

<10 ,6 g((o) 


< 10 12 g (2oo) 

< lO’ 16 g ((0) 


< 10' 16 g 


< lO' 10 g 


along spin axis 
within 1 urn 




Cryostat internal 
diameter > 1 m 



Centrifugal gradients due 
to skewed rotation axis 


< 10” g 


Centering of 
proof masses 
within 10 pm 


Use Niobium 
alloy blanket 
around detector. 
Degauss proof 


Proof masses Mass of cryostat 
with almost equal < 500 kg, 
moments of Intern, dia. > 1 m 
inertia. Belted 
cylinders not 
required. 


Centering of 
proof masses 
8* < 1 pm 



Limit magnetic 
moments outside 
sensor package 
M m < TBD A-m 2 
and r > TBD m 


N/A 



Verticality within 
0.1°; rate errors 
at release 
<0.1 °/s 


(0 = signal angular frequency; (flb = detector resonant angular frequency 


99 













































Experimental Activity (IFSI/CNR) 1 

Our partners have been working very hard in carrying out experimental activity on 
critical aspects of the detector development. The team at the Institute of Space Physics 
(IFSI/CNR) has already built a differential accelerometer prototype and carried out 
significant laboratory measurements on the prototype. 

The prototype was designed with the goal of exploring key aspects of a differential 
accelerometer and not for carrying out a preliminary test of the Equivalence Principle in 
the laboratory. As such, the prototype has two sensing masses of the same material with 
their centers of masses close together but not perfectly coincident. The prototype has 
several features in common with the instrument that we expect to develop for the flight 
experiment as follow: same capacitive pick-up system; same elastic suspension of the 
sensing masses; same measurement chain to extract the differential signal from the 
accelerometer. 

Figures 76 and 77 show a cross section of the differential accelerometer prototype 
when disassembled and assembled. Each sensing mass (blue and orange elements in the 
figure) is shaped to compenetrate (with leeway) into the opposing one so as to bring the 
two CMs close together (perfect coincidence is not a requisite for this prototype). 



Figure 6 Exploded view of differential accelerometer prototype 


Section contributed by V. Iafolla and S Nozzoli of IFSI/CNR funded through Italian Space Agency (ASI). 
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Figure 7 Cross section of (assembled) differential accelerometer prototype 



Figure 8 Simplified electrical diagram of signal pick-ups 

Each sensing mass can rotate (through elastic restraints) about the pivot axis on each 
side of the accelerometer. The two pivot axes are parallel to one another. The sensitive 
axis of the differential accelerometer is perpendicular to the pivot axes and to the 
longitudinal axis of the instrument (spin axis). Each sensing mass has two fixed 
capacitor plates for signal pickup. 


101 








MW 




tit im 




• ’> V • ‘jf* , 

. j ‘ > »n^V 




• <> *■' 

•• '^LfiT •" ' * “» 


>-•-» ■. • sn • v*v <• -L* ■£/ '> v v ,v ’ • ' * v '2L*o 

i > r *b-* ■ .s v - .. .: - . * - 


Figure 9 Pictures of differential accelerometer prototype 
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The output signals of the sensing masses are independently amplified and sent to the 
data acquisition system for conditioning, filtering and comparing. This prototype 
differential accelerometer builds on the heritage of (single) high-sensitivity 
accelerometers that were built by IFSI through the years. These accelerometers were 
tested extensively in the laboratory and in the field xn . Key characteristics of the 
differential accelerometer prototype are summarized in Table 1. 

Table 2 Key characteristics of differential accelerometer prototype 


Item 


Sensing mass 

0.22 kg 

Quality factor 

2900 

Resonance frequency 

18.12 Hz 

Preamplifier noise (K) 

0.76 

Preamplifier type 

AD743 (commercial) 



1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 

time (day) 


Figure 10 Tide signal filtered with 49-hour filter (1 mas = 2.78x1 0‘ 7 g) and measured 
with single accelerometer. August 1998. 
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The single accelerometer, in its various implementations, has demonstrated a 
remarkable sensitivity and repeatability. The accelerometer was able to resolve a 
modulated gravity signal (at 0.3 Hz) with a strength of less than 10' 9 g in the presence of 
seismic noise. Another single accelerometer is used to measure the Luni-Solar tides 
(working as a tilt meter) in die Gran Sasso laboratory (Italy). Excellent tidal signal over 
periods of months have been measured by the IFSI team. A month-long sample of a 
Luni-Solar tide in August 1998 is shown in Fig. 80 as an example of the performance of 
the single accelerometer in the field. 

Damping of transient oscillations 

The differential accelerometer prototype was developed and built to test, thus far, key 
aspects of differential acceleration measurements as follows: (a) abatement of the natural 
dynamics excited by the instrument release into the capsule; and (b) rejection of the 
common-mode signals. Point (a) above is critical for the success of our experimental 
scheme. For an instrument that is required to measure acceleration of 10' 14 g, the release 
is an abrupt event that pushes the detector well beyond its saturation point. Moreover, a 
(flight) instrument with a Quality factor (Q) of order 10 5 and a resonance frequency of a 
few Hz, it would take a very long time for the oscillation to decay to within the 
instrument’s dynamic range. The strategy that we plan to follow in the experiment is to 
reduce the Q factor during the first few seconds of fall from 10 5 down to a few units. The 
reduction of the Q factor is accomplished by inserting a resistance in the feedback control 
loops of the accelerometer. The technique has been tested successfully on the prototype 
accelerometer. 



Figure 11 Electrical diagram of one test mass pickup system with resistance added to the 
feedback loop 

Figure 82 shows the oscillation amplitude (from the oscilloscope output) of the 
accelerometer after an excitation. The accelerometer has a resonant frequency of 18.5 Hz 
and an (undamped) Q t = 2900. The introduction of 50x1 0 6 ohm resistive load in the 
feedback loop, reduces the Q from its undamped value to Q t ’ = 441 and, consequently, 
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the oscillation amplitude decays much more rapidly. Subsequently, the resistance is 
removed to demonstrate that this operation does not re-excite the oscillation amplitude. 


a X 6.75 V e. 134048 Tine 



Figure 12 Oscillation amplitude decay vs. time with and without resistance inserted in 
the feedback loop 

It is worth re mindin g that the decay of the proof mass oscillations is driven by the 
total quality factor of the electro-mechanical system which is as follows: 


Qt 



Qm Qe 


with 


1 o <*>qRC 
Qe 1 + (a) 0 RC) 2 


where P is the electro-mechanical coupling factor 




CE 2 



that is the ratio of electrical energy to mechanical energy. The electro-mechanical 
coupling factor of the instrument prototype is a low P = 0.01. The flight instrument will 
have a value much closer to unity thanks to a lower value of Wo and higher value of the 
capacitance C. Consequently, an electrical resistance of the value adopted for this test 
coupled into the feed-back loop of the flight instrument will imply a reduction of the 
value of Q t to a few units. In summary, this damping technique is able to provide the 
desired damping performance. 
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Common-mode rejection factor 

One important characteristic of a differential accelerometer is its ability to reject 
perturbations that are not differential, i.e., common-mode disturbances. This ability is 
quantified by the common-mode rejection factor (CMRF). 



Figure 13 Experimental setup to perturb differential accelerometer with periodic 
common-mode disturbances. 
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The IFSI laboratory has built a rotating system whose rotation axis can be tilted away 
from the vertical (see Fig. 83). The accelerometer is mounted inside the rotating device 
with the sensitive axis orthogonal to the spin axis. The rotating outfit is tilted by about 
10" 3 rad away from the vertical and spun at a slow rotation with a frequency of 0.15 Hz. 
As the differential accelerometer slowly rotates, it measures a small component 
(proportional to the tilt) of the Earth’s gravity field at the rotation frequency. This 
perturbation is a common-mode disturbance which affects equally the two proof masses. 




Figure 14 Accelerometer outputs: (a) single acceleration from proof mass 1 and 2 and 
(b) differential acceleration. 

The differential output generated by the calibration signal is then treated by the 
software through spectral analysis of the individual signals to adjust the proportionality 
factors of the two sensing masses and the phase differences of their responses. 
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* Figure 15 Spectra of individual and differential acceleration outputs: (a) after amplitude 

calibration only and (b) after amplitude and phase calibration. 

Figure 85 shows that after calibrating for amplitude and phase a 10 4 attenuation is 
readily obtained for the differential signal. This level of attenuation is effective not only 
at the perturbation frequency of 0.15 Hz but also over a larger frequency band. An 
attenuation of 10 4 or equivalently a common-mode rejection factor of 10"* meets the 
present requirement on the CMRF for the proposed tests of the Equivalence Principle. 
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Concluding remarks 

The laboratory activity consisted in the construction of a laboratory prototype of a 
differential accelerometer. The laboratory prototype has been used to conduct key tests 
on the differential instrument. We demonstrated the ability to damp quickly transient 
oscillations by utilizing a resistive load in the feedback loops and then removing that load 
to reestablish a high quality factor of the detector. A rotating divide with tilt control was 
also built. This device was utilized to impart (through the Earth’s gravity) common- 
mode perturbations to the differential accelerometer. These calibration disturbances have 
been used to trim the acceleration outputs of the individual proof masses in order to 
obtain a common-mode rejection factor better than 1 O' 4 in a sufficiently large frequency 
band centered at the spin frequency. 


\» 
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Key Areas to be Analyzed 


Instrument package 

Simulate dynamics of instrument package and differential accelerometer during 
free fall in the presence of the gravity gradient field. 

Extract from simulated data an EP violation signal at the sensitivity threshold in 
the presence of errors at release, CM position errors and intrinsic noise of the 
detector. 

Analyze effect of gravity quadrupole disturbance on instrument package and 
compare results with requirement on common-mode rejection factor. 

Select inertia characteristics of instrument package and derive balancing 
requirements. 

Leveling/Release mechanism 

Develop dynamic model of capsule attached to balloon gondola and study 
dynamics and stabilization before double release. Fine tune design concept. 

Cryostat design (in cooperation with Janis Research and Topsfield Engineering) 
Dewar layout. Define interfaces. Evaluate mass and cost of cryostat. 

Detector design 

Cooperate with our non-US partners on laboratory testing of instrument prototype 
to address critical issues related to the flight experiment; 

Cooperate with our non-US partners on proof mass design (for flight instrument) 
and mechanization of the detector. 

Magnetic disturbances 

Revisit the early estimates of magnetic disturbances as a result of updated 
configuration of instrument package. 

Derive requirements of magnetic cleanliness for typical levels of ferromagnetic 
impurities in the sensing masses. 
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Appendix A: Implementation and Computer Code 

The purpose of the code is to compute the force and the torque acting on each test 
mass, given arbitrary capsule. The test mass is characterized by its mass, and moment of 
inertia distribution (first, second and third order). The capsule is decomposed into few 
basic objects, such as cylinder, caps, and point mass, etc. The position of a test body * 

relative to the capsule is given by its center of mass, C.M., and by its rotation matrix, O . 

The radius vector of the center of mass is expressed in capsule frame. The rotation matrix k 

defines the rotation between the capsule frame, and the body frame (not necessarily the 
principal frame). The origin of the capsule coordinate system is at the center of the 
cylinder base, where the Z-axis (axisymmetry axis) is pointing upward. 

The code was designed as general and as flexible as possible, and also with 
computational efficiency in mind. The user may define and add any number of basic 
capsule objects. The code was written in Matlab, and it contains m files as source code 
and mat files as data files. It also contains script language (in Matlab), for the purpose of 
automatic generation and reading of files. 

The following figure shows a schematic description of the code. Attached below is 
the Matlab help file, located at the code directory. It describes the code and gives specific 
instructions how to run the code. 
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G_help.txt 


General description 


The distributed mass elements for each shape are defined (by user) in 
the files: shape#. m 

Each file defines a particular shape (cylinder, cup, point mass, etc) . 
The output data is the element 

coordinates vector: x#, y#, z#, and a mass vector: m# . These are saved 
into shape#. mat files. 

Note that the shape enumeration does not have to be in a specific 
order. The user may add any number of 

shape files, and these should save the date to the appropriate 
shape#. mat files. 

The test body is defined in the file: test_mass# .m, where # is an index 
(ex: 1 or 2 ) . 

The user defined parameters are the body index, the mass and the 
principal moments of inertia of the body. 

The results are saved into the files: test_mass# .mat 

The file: parameters. m defined (by user) global parameters. 

It enables the user to include any combination of a given shape from 
shape#. dat by defining the 

index array: i_shape. This array includes the shape number, (ex: 
i_shape= [1 3 6 ] ) . 

The file lumped_mass . m reads any shape#. mat specified by i_shape, and 
aggravates the coordinates vectors 

and the mass vectors into single vectors: x,y,z and m. These vectors 
are saved into the file: shapes. mat 

The force/torque computation is done in the function file: G_comp.m 
It loads the data files: shapes. mat and test_mass . mat 

The input arguments to this function are: 1. body. It is a cell array 
contains the CM position and the 

transformation matrix for the orientation of the body, relative to the 
capsule frame. 2. The body index. 

The output of this function are arrays of force and array of torque. 
Software usage 


The package include the files: 

initialize .m 
parameters. m 
test^raass .m 
shape# . m .... 

* G_comp . m 

G_driver .m 

These files generate a stand alone application, when using the driver: 
G_driver to run the code. 

However, for a future usage, the function: G_comp can be a part of a 
dynamical simulation code, 

with the interface of the input/output of this function. 

The software consists of two functionalities: 
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1. Initialization of the data files: *.mat and the global 
parameters . 

2. Gravitational computation. 

Once the user initializes, G_comp.m can be run independently. 

To run the code, the following action should be taken: 

- Defining the parameters in: shape#. m ..., test__mass.m and in 
parameters . m 

- Typing as command line: initialize 

- Defining in G_driver: body index, body{l}=CM and body { 2 } =Trans , or 

defining it in other calling routine. * 

- Typing as command line: [ F, T] =G_comp (body, index) 


i' 




i 
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Appendix B: Gravitational Model for the Caps 


The discretization of the flat cup is by using a rectangular grid, where the lumped- 
masses are allocated at the mid-layer. 

In order to verify the precision of the discrete model, we will compare it with an 
analytical model. 

It is easy to see that the axial force between a point mass, HI , and a finite size disk 
with radius R c and density p , is: 

Z+t 1 2 Rc rii 

Fz = 2jiGpm f f — dr dll 

ij=z-t/2r=olr~ + rj'J 

where Z is the radius vector between the point mass and the center of the disk. 

The resulting force is: 


Fz 


= 2 ji Gpm 
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However, our model is equivalent with a thin layer distribution of masses. So, the 
analytical expression for the force between a point mass and a thin disk is computed by 
applying the limiting process for t : 


Fz = 2jx Gm lim 
t^O r->0 t 



1 1 V 

t~i 

Z+-t\ 

V 

\ 2 ) 



where p ^ = t ’ p , thus p ^ is the area density. The resulting approximation is: 


/ 

Fz = 2 JtGpmt 1 

t-* 0 

V 


V 


z2+R c , 


The question now is where is the best location of the thin disk (i.e. somewhere 
between the lower and upper surfaces of the cup). Figure Cl compares the axial force due 
to a finite thickness disk with the force due to a thin disk. The location of the thin disk 
varies between the lower surface and the upper surface. It is normalized by the thickness, 
so its value at the lower surface is 0 and its value at the upper surface is 1 . Z is the radius 
vector between a point mass of 1 kg, and the middle of the cup. It is obvious that the best 
location for the thin disk is at the middle of the cup. 
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This argument can be proved analytically as well. Suppose that Z q is the radius 

vector for the thin disk (unknown, so far), and Z is the radius vector to the middle of the 
cup. Comparing the corresponding forces results in, 



This equation can be solved for t/R c « 1, and a straightforward expansion leads to. 



Normalized Cup Thickness 


Figure B 1 . Forces due to various locations of thin disk compared with forces due to 
finite thickness disk 

The next question concerns the sufficient discretization. Let N be the number of grids 
along the diameter. (The number of cells, that is the number of lumped-masses, becomes 

ji 2 

closer to — /V as N~- * 0 ° ). Figure C2 shows the axial force between the two caps, for a 

4 
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point mass of 1 kg. Even a crude discretization (N=4) gives a reasonable approximation. 
In our code we take N=20; sufficiently to represent the caps. 



Z[m] 

Figure B2. Effect of discretization on axial force between two caps 
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Appendix C: Gravitational Model for the Cylinder 


The discretization of the flat cup is by done by dividing the circumference into 
sectors, and the height into rectangular grid, in a way that the horizontal and the vertical 
solid angles are equal. 

In order to verify the precision of the discrete model, we will compare it with an 
analytical model. 

The axial force between a point mass, m, at height h from the cylinder button, and a 
finite thickness cylinder with internal radius /?,, external radius R 2 , height H and 
density p is: 

H-h Ri 

Fz = 2 jt Gp m f f 

il=-li r=R ] 

The first integration is related to the force between the mass element and a hollow 
disk, at a distance TJ . The second integration is along the height of the cylinder. 

The resulting force is: 

Fz = 2jtGpm[jRf+(H-h) 2 ~^R 2 +h 2 -^R^HH-h) 2 + ^R 2 +h 2 ] 

Note that in order to work with coordinate system that is centered at the middle height 
( X,Y,Z), we have to apply the following transformation: h*—H / 2+Z 

Next, we compute the thin cylinder approximation force (as our discrete model), and 
compare it with the finite thickness force. 


Fz = 2 jzGm lim ^-[ Jr? +(H-h ) 2 - J^+h 2 - ^(R]+t) 2 +{H -h) 2 - + y\(R\+t) 2 +h 2 
0 t -*0 1 v v 



The resulting approximation is: 


/ 


Fz 

t-* 0 


= 2 JiGmpt 



1 
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The best R is found by comparing the thin cylinder approximation with the finite 
thickness formula. Figure Cl suggests that R = (/?, + R 2 )/2 is the best place to locate the 
approximate cylinder. 


*1 
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Normalized Cylinder Thickness 


Figure Cl : Axial force due to various radius of thin cylinder compared with 
an axial force due to finite thickness cylinder 

The accuracy of the discretization is demonstrated in figure C2. The axial force 
computed by the thin cylinder approximation is compared with numerical solution of the 
discrete model. The angular interval between the lumped-masses is 1 degree. Note that 
even an angular interval of 6 degrees gives a reasonable model (the maximal error is less 
than 1%). 
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Figure C2. Axial force for a test mass along the cylinder axisymmetrical axis: analytical 
Vs numerical solutions 




Appendix D: Interpolation of the Gravitational Force Field by Radial Basis 
Functions 


* 


I 
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The computation of the gravitational force inside the capsule is time consuming, 
especially when considering a dynamical simulation of the test masses. In this case, the 
gravitational field has to be evaluated at each integration instance. The computational 
efficiency may increases dramatically if we represent the gravitational field by a low 
order functional approximation. 

Our interpolation method is generally subjected to two requirements. First, the 
approximate surface should be close enough to the computational surface. Second, the 
approximate surface should be as smooth as possible. The second requirement guarantees 
the computation of a gravity gradient field. Therefore, we avoid using any interpolation 
that introduces high frequencies and has ill convergence (for example, polynomial 
approximation). A feasible interpolation method is the following Radial Basis Functions 
(RBF) method. 

Consider a mapping from a d-dimensional input space, X, to a one-dimensional target 
space, T. In our problem the input space consists of the radial direction, x, and the axial 
direction, z, (d = 2). The target space is a component of the gravitational force, Fx or Fz. 
The data set consists of N input vectors X, together with the corresponding N targets. The 
goal is to find a function S(X), such that S(X n ) = t n ; n = 1 ,..., N . That is, the mapping 
should be exact at the fitting points, and approximates the rest of the field. 

The RBF approach introduces a set of N basis functions, one for each data point. We 
adopted RBF in the Gaussian form: 




= exp 


2 a; 


Note that each training data serves as a center for the basic function, R n . 

The variance o determines the locality of each RBF. Small variance will result is a 
good local fitting but worse fitting between the training data, and worse smoothing. 
While a large variance will result in a better smoothing but worse local fitting. The 
interpolation function is a linear combination of the RBF, in the form: 


S(X) = Jw exp 


(*-U 2 

2 a; 


Assuming constant variance, the mapping relations of the N fitting points results in 
the following linear system: 
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The matrix $ is well conditioned under certain conditions, determine by the size of 
the system and the covariance. The solution for the weights is: 

w = <r' T . 

The RBF for the force field was carried out as follows. The training data is taken in 
the region of the expected trajectory (21 points), as demonstrated in the figure below. 
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Figure Dl. The training data for the RBF approximation inside a cylinder 
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The forces Fx and Fz were computed for each training point. The quality of the 
fitness is demonstrated for Fx (harder surface to fit than Fz). The variance was taken as 
lm. The worst fitting is about 0.2%. 


xIO 9 



- 1.5 -1 - 0.5 0 0.5 1 1.5 

Z[m] 


Figure D2. Validation of the RBF approximation for the radial force 

The interpolated surfaces for Fx and Fz are shown below. The computation of these is 
very efficient now, because it requires the evaluation of only 21 exponential terms. Note 
that this example is related to a simple open cylinder. The force field due to a more 
complicated chamber, with additional attracting elements, will be handled with almost 
the same computational effort. That is, after the computation of the force field for the 
training points, only 21 RBF evaluations are required for approximating the force field at 
a different location. 
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Figure D3: RBF approximate surfaces for Fx and Fz 


Note that there is a fitting problem when the force field consists of terms of different 
order of magnitude. In this case the RBF will fit only the largest terms. Given a force 

field expressed as: F = F ^ +• • •, a feasible solution is to approximate each term 

by a different RBF, such that the approximate force field becomes: 

F = F (1) + F& +.... 

Another approach is to start from the gravitational potential. After evaluating the 
potential of the capsule/test-body at the training points, and computing the RBF weights, 
taking the gradient of the RBF representation gives the forces. A similar approach, for 
evaluating the gravity gradient, is demonstrated in the following section. 


v 


One of the advantages of using RBF is the smoothness that enables us to compute the 
gradient. The gravity gradient field can be also interpolated by analytical computation, 
and then approximating it by RBF. However, the analytical expression for the gravity 
gradient of a high-order potential may be a very complicated. Therefore, our approach is 
to compute the RBF for the force field and use it as a generating function for the gravity 
gradient. 
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Let F be the force field expanded in RBF. 
F = ^ W n exp ,J — 


where r = (x-x„) : +(y-y„) 2 + (z-z„) 2 , 
F = [Fx Fy Fz] T . 

The gravity gradient is then, 

v x F x v y f x v z f x ' 

G = S7 r F= V X Fy VyFy V z F Y 
v X F Z VyF Z ^Z F Z 


W , and 


Given a RBF force field at a point r, the RBF force field at an adjacent point is: 
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Let cp(r n ) = exp — be a local radial basis function, then its gradient is: 
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The gravity gradient expressed in RBF is the following superposition of the gradients 
of the local radial basis functions. 

G-ilVy r cp(r,) 

M«l 

The final form of the RBF gravity gradient is: 

w'(x-x ) w ' ( v - V ) vv'fz-z ) 

1 v r~ 

G = — 7 2 w '.( x ~ x » ) O' - y . ) w ,! ( z - z „ ) ex P - zr^ 
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